Replace sysimage.jl by juliaw; add deps/formatter

Improve benchmark scripts
Improve sysimage.jl
2025-12-06 08:18:51 -06:00 · 2021-08-31 09:51:36 -05:00 · 2021-08-31 08:03:21 -05:00 · 2021-08-20 04:51:09 -05:00 · 2021-08-19 07:07:25 -05:00 · 2021-08-05 17:04:37 -05:00
49 changed files with 856 additions and 1355 deletions
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,4 @@ TODO.md
 docs/_build
 .vscode
 Manifest.toml
+*/Manifest.toml
--- a/23
+++ b/23
@@ -2,31 +2,22 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-JULIA := julia --color=yes --project=@.
 VERSION := 0.2

-build/sysimage.so: src/utils/sysimage.jl Project.toml Manifest.toml
-	mkdir -p build
-	mkdir -p benchmark/results/test
-	cd benchmark; $(JULIA) --trace-compile=../build/precompile.jl benchmark.jl test/case14
-	$(JULIA) src/utils/sysimage.jl
-
 clean:
-	rm -rf build/*
+	rm -rfv build

 docs:
 	cd docs; make clean; make dirhtml
 	rsync -avP --delete-after docs/_build/dirhtml/ ../docs/$(VERSION)/
 	
-test: build/sysimage.so
-	@echo Running tests...
-	$(JULIA) --sysimage build/sysimage.so -e 'using Pkg; Pkg.test("UnitCommitment")' | tee build/test.log
-
-
 format:
-	julia -e 'using JuliaFormatter; format(["src", "test", "benchmark"], verbose=true);'
+	cd deps/formatter; ../../juliaw format.jl

-install-deps:
-	julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.14.4"))'
+test: test/Manifest.toml
+	./juliaw test/runtests.jl
+
+test/Manifest.toml: test/Project.toml
+	julia --project=test -e "using Pkg; Pkg.instantiate()"

 .PHONY: docs test format install-deps
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.2.2"

 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
@@ -15,10 +16,10 @@ 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]
-Cbc = "0.7"
 DataStructures = "0.18"
 Distributions = "0.25"
 GZip = "0.5"
@@ -27,11 +28,3 @@ JuMP = "0.21"
 MathOptInterface = "0.9"
 PackageCompiler = "1"
 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", "Gurobi"]
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+DocOpt = "968ba79b-81e4-546f-ab3a-2eecfa62a9db"
 Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
--- a/benchmark/benchmark.jl
+++ b/benchmark/benchmark.jl
@@ -1,158 +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 Distributed
-using Pkg
-Pkg.activate(".")
-
-@everywhere using Pkg
-@everywhere Pkg.activate(".")
-
-@everywhere using UnitCommitment
-@everywhere using JuMP
-@everywhere using Gurobi
-@everywhere using JSON
-@everywhere using Logging
-@everywhere using Printf
-@everywhere using LinearAlgebra
-@everywhere using Random
-
-@everywhere import UnitCommitment:
-    ArrCon2000,
-    CarArr2006,
-    DamKucRajAta2016,
-    Formulation,
-    Gar1962,
-    KnuOstWat2018,
-    MorLatRam2013,
-    PanGua2016,
-    XavQiuWanThi2019
-
-@everywhere UnitCommitment._setup_logger()
-
-function main()
-    cases = [
-        "pglib-uc/ca/2014-09-01_reserves_0",
-        "pglib-uc/ca/2014-09-01_reserves_1",
-        "pglib-uc/ca/2015-03-01_reserves_0",
-        "pglib-uc/ca/2015-06-01_reserves_0",
-        "pglib-uc/ca/Scenario400_reserves_1",
-        "pglib-uc/ferc/2015-01-01_lw",
-        "pglib-uc/ferc/2015-05-01_lw",
-        "pglib-uc/ferc/2015-07-01_hw",
-        "pglib-uc/ferc/2015-10-01_lw",
-        "pglib-uc/ferc/2015-12-01_lw",
-        "pglib-uc/rts_gmlc/2020-04-03",
-        "pglib-uc/rts_gmlc/2020-09-20",
-        "pglib-uc/rts_gmlc/2020-10-27",
-        "pglib-uc/rts_gmlc/2020-11-25",
-        "pglib-uc/rts_gmlc/2020-12-23",
-        "or-lib/20_0_1_w",
-        "or-lib/20_0_5_w",
-        "or-lib/50_0_2_w",
-        "or-lib/75_0_2_w",
-        "or-lib/100_0_1_w",
-        "or-lib/100_0_4_w",
-        "or-lib/100_0_5_w",
-        "or-lib/200_0_3_w",
-        "or-lib/200_0_7_w",
-        "or-lib/200_0_9_w",
-        "tejada19/UC_24h_290g",
-        "tejada19/UC_24h_623g",
-        "tejada19/UC_24h_959g",
-        "tejada19/UC_24h_1577g",
-        "tejada19/UC_24h_1888g",
-        "tejada19/UC_168h_72g",
-        "tejada19/UC_168h_86g",
-        "tejada19/UC_168h_130g",
-        "tejada19/UC_168h_131g",
-        "tejada19/UC_168h_199g",
-    ]
-    formulations = Dict(
-        "Default" => Formulation(),
-        "ArrCon2000" => Formulation(ramping = ArrCon2000.Ramping()),
-        "CarArr2006" => Formulation(pwl_costs = CarArr2006.PwlCosts()),
-        "DamKucRajAta2016" =>
-            Formulation(ramping = DamKucRajAta2016.Ramping()),
-        "Gar1962" => Formulation(pwl_costs = Gar1962.PwlCosts()),
-        "KnuOstWat2018" =>
-            Formulation(pwl_costs = KnuOstWat2018.PwlCosts()),
-        "MorLatRam2013" => Formulation(ramping = MorLatRam2013.Ramping()),
-        "PanGua2016" => Formulation(ramping = PanGua2016.Ramping()),
-    )
-    trials = [i for i in 1:5]
-    combinations = [
-        (c, f.first, f.second, t) for c in cases for f in formulations for
-        t in trials
-    ]
-    shuffle!(combinations)
-    @sync @distributed for c in combinations
-        _run_combination(c...)
-    end
-end
-
-@everywhere function _run_combination(
-    case,
-    formulation_name,
-    formulation,
-    trial,
-)
-    name = "$formulation_name/$case"
-    dirname = "results/$name"
-    mkpath(dirname)
-    if isfile("$dirname/$trial.json")
-        @info @sprintf("%-4s %-16s %s", "skip", formulation_name, case)
-        return
-    end
-    @info @sprintf("%-4s %-16s %s", "run", formulation_name, case)
-    open("$dirname/$trial.log", "w") do file
-        redirect_stdout(file) do
-            redirect_stderr(file) do
-                return _run_sample(case, formulation, "$dirname/$trial")
-            end
-        end
-    end
-    @info @sprintf("%-4s %-16s %s", "done", formulation_name, case)
-end
-
-@everywhere function _run_sample(case, formulation, prefix)
-    total_time = @elapsed begin
-        @info "Reading: $case"
-        time_read = @elapsed begin
-            instance = UnitCommitment.read_benchmark(case)
-        end
-        @info @sprintf("Read problem in %.2f seconds", time_read)
-        BLAS.set_num_threads(4)
-        model = UnitCommitment.build_model(
-            instance = instance,
-            formulation = formulation,
-            optimizer = optimizer_with_attributes(
-                Gurobi.Optimizer,
-                "Threads" => 4,
-                "Seed" => rand(1:1000),
-            ),
-            variable_names = true,
-        )
-        @info "Optimizing..."
-        BLAS.set_num_threads(1)
-        UnitCommitment.optimize!(
-            model,
-            XavQiuWanThi2019.Method(time_limit = 3600.0, gap_limit = 1e-4),
-        )
-    end
-    @info @sprintf("Total time was %.2f seconds", total_time)
-    @info "Writing solution: $prefix.json"
-    solution = UnitCommitment.solution(model)
-    UnitCommitment.write("$prefix.json", solution)
-    @info "Verifying solution..."
-    return UnitCommitment.validate(instance, solution)
-    # @info "Exporting model..."
-    # return JuMP.write_to_file(model, model_filename)
-end
-
-if length(ARGS) > 0
-    _run_sample(ARGS[1], UnitCommitment.Formulation(), "tmp")
-else
-    main()
-end
--- a/benchmark/run.jl
+++ b/benchmark/run.jl
@@ -0,0 +1,207 @@
+# 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.
+
+doc = """UnitCommitment.jl Benchmark Runner
+
+Usage:
+  run.jl [-s ARG]... [-m ARG]... [-c ARG]... [-f ARG]... [options]
+
+Examples:
+
+    1. Benchmark all solvers, methods and formulations:
+
+        julia run.jl
+
+    2. Benchmark formulations "default" and "ArrCon200" using Gurobi:
+
+        julia run.jl -s gurobi -f default -f ArrCon2000
+
+    3. Benchmark a few test cases, using all solvers, methods and formulations:
+
+        julia run.jl -c or-lib/20_0_1_w -c matpower/case1888rte/2017-02-01
+
+    4. Solve 4 test cases in parallel, with 2 threads available per worker:
+
+        JULIA_NUM_THREADS=2 julia --procs 4 run.jl
+
+Options:
+  -h --help             Show this screen.
+  -s --solver=ARG       Mixed-integer linear solver (e.g. gurobi)
+  -c --case=ARG         Unit commitment test case (e.g. or-lib/20_0_1_w)
+  -m --method=ARG       Solution method (e.g. default)
+  -f --formulation=ARG  Formulation (e.g. ArrCon2000)
+  --time-limit=ARG      Time limit in seconds [default: 3600]
+  --gap=ARG             Relative MIP gap tolerance [default: 0.001]
+  --trials=ARG          Number of trials [default: 5]
+"""
+
+using Distributed
+using Pkg
+Pkg.activate(".")
+@everywhere using Pkg
+@everywhere Pkg.activate(".")
+
+using DocOpt
+args = docopt(doc)
+
+@everywhere using UnitCommitment
+@everywhere UnitCommitment._setup_logger()
+
+using UnitCommitment
+using Gurobi
+using Logging
+using JuMP
+
+import UnitCommitment:
+    ArrCon2000,
+    CarArr2006,
+    DamKucRajAta2016,
+    Formulation,
+    Gar1962,
+    KnuOstWat2018,
+    MorLatRam2013,
+    PanGua2016,
+    XavQiuWanThi2019
+
+# Benchmark test cases
+# -----------------------------------------------------------------------------
+cases = [
+    "pglib-uc/ca/2014-09-01_reserves_0",
+    "pglib-uc/ca/2014-09-01_reserves_1",
+    "pglib-uc/ca/2015-03-01_reserves_0",
+    "pglib-uc/ca/2015-06-01_reserves_0",
+    "pglib-uc/ca/Scenario400_reserves_1",
+    "pglib-uc/ferc/2015-01-01_lw",
+    "pglib-uc/ferc/2015-05-01_lw",
+    "pglib-uc/ferc/2015-07-01_hw",
+    "pglib-uc/ferc/2015-10-01_lw",
+    "pglib-uc/ferc/2015-12-01_lw",
+    "pglib-uc/rts_gmlc/2020-04-03",
+    "pglib-uc/rts_gmlc/2020-09-20",
+    "pglib-uc/rts_gmlc/2020-10-27",
+    "pglib-uc/rts_gmlc/2020-11-25",
+    "pglib-uc/rts_gmlc/2020-12-23",
+    "or-lib/20_0_1_w",
+    "or-lib/20_0_5_w",
+    "or-lib/50_0_2_w",
+    "or-lib/75_0_2_w",
+    "or-lib/100_0_1_w",
+    "or-lib/100_0_4_w",
+    "or-lib/100_0_5_w",
+    "or-lib/200_0_3_w",
+    "or-lib/200_0_7_w",
+    "or-lib/200_0_9_w",
+    "tejada19/UC_24h_290g",
+    "tejada19/UC_24h_623g",
+    "tejada19/UC_24h_959g",
+    "tejada19/UC_24h_1577g",
+    "tejada19/UC_24h_1888g",
+    "tejada19/UC_168h_72g",
+    "tejada19/UC_168h_86g",
+    "tejada19/UC_168h_130g",
+    "tejada19/UC_168h_131g",
+    "tejada19/UC_168h_199g",
+    "matpower/case1888rte/2017-02-01",
+    "matpower/case1951rte/2017-02-01",
+    "matpower/case2848rte/2017-02-01",
+    "matpower/case3012wp/2017-02-01",
+    "matpower/case3375wp/2017-02-01",
+    "matpower/case6468rte/2017-02-01",
+    "matpower/case6515rte/2017-02-01",
+]
+
+# Formulations
+# -----------------------------------------------------------------------------
+formulations = Dict(
+    "default" => Formulation(),
+    "ArrCon2000" => Formulation(ramping = ArrCon2000.Ramping()),
+    "CarArr2006" => Formulation(pwl_costs = CarArr2006.PwlCosts()),
+    "DamKucRajAta2016" => Formulation(ramping = DamKucRajAta2016.Ramping()),
+    "Gar1962" => Formulation(pwl_costs = Gar1962.PwlCosts()),
+    "KnuOstWat2018" => Formulation(pwl_costs = KnuOstWat2018.PwlCosts()),
+    "MorLatRam2013" => Formulation(ramping = MorLatRam2013.Ramping()),
+    "PanGua2016" => Formulation(ramping = PanGua2016.Ramping()),
+)
+
+# Solution methods
+# -----------------------------------------------------------------------------
+methods = Dict(
+    "default" =>
+        XavQiuWanThi2019.Method(time_limit = 3600.0, gap_limit = 1e-4),
+)
+
+# MIP solvers
+# -----------------------------------------------------------------------------
+optimizers = Dict(
+    "gurobi" => optimizer_with_attributes(
+        Gurobi.Optimizer,
+        "Threads" => Threads.nthreads(),
+    ),
+)
+
+# Parse command line arguments
+# -----------------------------------------------------------------------------
+if !isempty(args["--case"])
+    cases = args["--case"]
+end
+if !isempty(args["--formulation"])
+    formulations = filter(p -> p.first in args["--formulation"], formulations)
+end
+if !isempty(args["--method"])
+    methods = filter(p -> p.first in args["--method"], methods)
+end
+if !isempty(args["--solver"])
+    optimizers = filter(p -> p.first in args["--solver"], optimizers)
+end
+const time_limit = parse(Float64, args["--time-limit"])
+const gap_limit = parse(Float64, args["--gap"])
+const ntrials = parse(Int, args["--trials"])
+
+# Print benchmark settings
+# -----------------------------------------------------------------------------
+function printlist(d::Dict)
+    for key in keys(d)
+        @info "  - $key"
+    end
+end
+
+function printlist(d::Vector)
+    for key in d
+        @info "  - $key"
+    end
+end
+
+@info "Computational environment:"
+@info "  - CPU: $(Sys.cpu_info()[1].model)"
+@info "  - Logical CPU cores: $(length(Sys.cpu_info()))"
+@info "  - System memory: $(round(Sys.total_memory() / 2^30, digits=2)) GiB"
+@info "  - Available workers: $(nworkers())"
+@info "  - Available threads per worker: $(Threads.nthreads())"
+
+@info "Parameters:"
+@info "  - Number of trials: $ntrials"
+@info "  - Time limit (s): $time_limit"
+@info "  - Relative MIP gap tolerance: $gap_limit"
+
+@info "Solvers:"
+printlist(optimizers)
+
+@info "Methods:"
+printlist(methods)
+
+@info "Formulations:"
+printlist(formulations)
+
+@info "Cases:"
+printlist(cases)
+
+# Run benchmarks
+# -----------------------------------------------------------------------------
+UnitCommitment._run_benchmarks(
+    cases = cases,
+    formulations = formulations,
+    methods = methods,
+    optimizers = optimizers,
+    trials = 1:ntrials,
+)
--- a/deps/formatter/Project.toml
+++ b/deps/formatter/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
+
+[compat]
+JuliaFormatter = "0.14.4"
--- a/deps/formatter/format.jl
+++ b/deps/formatter/format.jl
@@ -0,0 +1,9 @@
+using JuliaFormatter
+format(
+    [
+        "../../src", 
+        "../../test",
+        "../../benchmark/run.jl",
+    ],
+    verbose=true,
+)
--- a/68
+++ b/68
@@ -0,0 +1,68 @@
+#!/bin/bash
+# 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.
+
+if [ ! -e Project.toml ]; then
+    echo "juliaw: Project.toml not found"
+    exit 1
+fi
+
+if [ ! -e Manifest.toml ]; then
+    julia --project=. -e 'using Pkg; Pkg.instantiate()' || exit 1
+fi
+
+if [ ! -e build/sysimage.so -o Project.toml -nt build/sysimage.so ]; then
+    echo "juliaw: rebuilding system image..."
+
+    # Generate temporary project folder
+    rm -rf $HOME/.juliaw
+    mkdir -p $HOME/.juliaw/src
+    cp Project.toml Manifest.toml $HOME/.juliaw
+    NAME=$(julia -e 'using TOML; toml = TOML.parsefile("Project.toml"); "name" in keys(toml) && print(toml["name"])')
+    if [ ! -z $NAME ]; then
+        cat > $HOME/.juliaw/src/$NAME.jl << EOF
+module $NAME
+end
+EOF
+    fi
+
+    # Add PackageCompiler dependencies to temporary project
+    julia --project=$HOME/.juliaw -e 'using Pkg; Pkg.add(["PackageCompiler", "TOML", "Logging"])'
+
+    # Generate system image scripts
+    cat > $HOME/.juliaw/sysimage.jl << EOF
+using PackageCompiler
+using TOML
+using Logging
+
+Logging.disable_logging(Logging.Info)
+mkpath("$PWD/build")
+
+println("juliaw: generating precompilation statements...")
+run(\`julia --project="$PWD" --trace-compile="$PWD"/build/precompile.jl \$(ARGS)\`)
+
+println("juliaw: finding dependencies...")
+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])
+        println("  - \$(dep) [skip]")
+    else
+        println("  - \$(dep)")
+        push!(deps, Symbol(dep))
+    end
+end
+
+println("juliaw: building system image...")
+create_sysimage(
+    deps,
+    precompile_statements_file = "$PWD/build/precompile.jl",
+    sysimage_path = "$PWD/build/sysimage.so",
+)
+EOF
+    julia --project=$HOME/.juliaw $HOME/.juliaw/sysimage.jl $* 
+else
+    julia --project=. --sysimage build/sysimage.so $*
+fi
--- a/src/UnitCommitment.jl
+++ b/src/UnitCommitment.jl
@@ -48,8 +48,9 @@ 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("utils/benchmark.jl")
 include("validation/repair.jl")
 include("validation/validate.jl")

--- a/src/instance/read.jl
+++ b/src/instance/read.jl
@@ -266,15 +266,20 @@ function _from_json(json; repair = true)
    end

    instance = UnitCommitmentInstance(
-        T,
-        power_balance_penalty,
-        shortfall_penalty,
-        units,
-        buses,
-        lines,
-        reserves,
-        contingencies,
-        loads,
+        buses_by_name = Dict(b.name => b for b in buses),
+        buses = buses,
+        contingencies_by_name = Dict(c.name => c for c in contingencies),
+        contingencies = contingencies,
+        lines_by_name = Dict(l.name => l for l in lines),
+        lines = lines,
+        power_balance_penalty = power_balance_penalty,
+        price_sensitive_loads_by_name = Dict(ps.name => ps for ps in loads),
+        price_sensitive_loads = loads,
+        reserves = reserves,
+        shortfall_penalty = shortfall_penalty,
+        time = T,
+        units_by_name = Dict(g.name => g for g in units),
+        units = units,
    )
    if repair
        UnitCommitment.repair!(instance)
--- a/src/instance/structs.jl
+++ b/src/instance/structs.jl
@@ -69,17 +69,21 @@ mutable struct PriceSensitiveLoad
    revenue::Vector{Float64}
 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}
+Base.@kwdef mutable struct UnitCommitmentInstance
+    buses_by_name::Dict{AbstractString,Bus}
    buses::Vector{Bus}
-    lines::Vector{TransmissionLine}
-    reserves::Reserves
+    contingencies_by_name::Dict{AbstractString,Contingency}
    contingencies::Vector{Contingency}
+    lines_by_name::Dict{AbstractString,TransmissionLine}
+    lines::Vector{TransmissionLine}
+    power_balance_penalty::Vector{Float64}
+    price_sensitive_loads_by_name::Dict{AbstractString,PriceSensitiveLoad}
    price_sensitive_loads::Vector{PriceSensitiveLoad}
+    reserves::Reserves
+    shortfall_penalty::Vector{Float64}
+    time::Int
+    units_by_name::Dict{AbstractString,Unit}
+    units::Vector{Unit}
 end

 function Base.show(io::IO, instance::UnitCommitmentInstance)
--- a/src/model/build.jl
+++ b/src/model/build.jl
@@ -21,8 +21,6 @@ Arguments
 - `optimizer`:
    the optimizer factory that should be attached to this model (e.g. Cbc.Optimizer).
    If not provided, no optimizer will be attached.
- `formulation`:
-    the details of which constraints, variables, etc. to use
 - `variable_names`: 
    If true, set variable and constraint names. Important if the model is going
    to be exported to an MPS file. For large models, this can take significant
--- a/src/model/formulations/ArrCon2000/ramp.jl
+++ b/src/model/formulations/ArrCon2000/ramp.jl
@@ -2,15 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_ramp_eqs!(model, unit, formulation_prod_vars, formulation_ramping, formulation_status_vars)::Nothing
-
-Ensure constraints on ramping are met.
-Based on Arroyo and Conejo (2000).
-Eqns. (24), (25) in Knueven et al. (2020).
-
-Adds constraints identified by `ArrCon200.Ramping` to `model` using variables `Gar1962.ProdVars` and `is_on` from `Gar1962.StatusVars`.
-"""
 function _add_ramp_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -31,6 +22,7 @@ function _add_ramp_eqs!(
    reserve = model[:reserve]
    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]
@@ -40,8 +32,6 @@ function _add_ramp_eqs!(
    switch_off = model[:switch_off]
    switch_on = model[:switch_on]

-    is_initially_on = _is_initially_on(g) > 0
-
    for t in 1:model[:instance].time
        # Ramp up limit
        if t == 1
@@ -66,7 +56,7 @@ function _add_ramp_eqs!(
            min_prod_last_period =
                g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]

-            # Equation (24) in Knueven et al. (2020)
+            # Equation (24) in Kneuven et al. (2020)
            eq_ramp_up[gn, t] = @constraint(
                model,
                max_prod_this_period - min_prod_last_period <=
@@ -97,7 +87,7 @@ function _add_ramp_eqs!(
            min_prod_this_period =
                g.min_power[t] * is_on[gn, t] + prod_above[gn, t]

-            # Equation (25) in Knueven et al. (2020)
+            # Equation (25) in Kneuven et al. (2020)
            eq_ramp_down[gn, t] = @constraint(
                model,
                max_prod_last_period - min_prod_this_period <=
--- a/src/model/formulations/ArrCon2000/structs.jl
+++ b/src/model/formulations/ArrCon2000/structs.jl
@@ -13,12 +13,6 @@ module ArrCon2000

 import ..RampingFormulation

-"""
-Constraints
---
-* `eq_ramp_up`: Equation (24) in Knueven et al. (2020)
-* `eq_ramp_down`: Equation (25) in Knueven et al. (2020)
-"""
 struct Ramping <: RampingFormulation end

 end
--- a/src/model/formulations/CarArr2006/pwlcosts.jl
+++ b/src/model/formulations/CarArr2006/pwlcosts.jl
@@ -2,12 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_production_piecewise_linear_eqs!
-
-Ensure respect of production limits along each segment.
-Creates constraints `CarArr2006.PwlCosts` using variables `Gar1962.StatusVars`
-"""
 function _add_production_piecewise_linear_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -27,7 +21,7 @@ function _add_production_piecewise_linear_eqs!(
    for t in 1:model[:instance].time
        gn = g.name
        for k in 1:K
-            # Equation (45) in Knueven et al. (2020)
+            # Equation (45) in Kneuven et al. (2020)
            # NB: when reading instance, UnitCommitment.jl already calculates
            #     difference between max power for segments k and k-1 so the
            #     value of cost_segments[k].mw[t] is the max production *for
@@ -42,14 +36,14 @@ function _add_production_piecewise_linear_eqs!(
            set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])

            # Definition of production
-            # Equation (43) in Knueven et al. (2020)
+            # Equation (43) in Kneuven et al. (2020)
            eq_prod_above_def[gn, t] = @constraint(
                model,
                prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
            )

            # Objective function
-            # Equation (44) in Knueven et al. (2020)
+            # Equation (44) in Kneuven et al. (2020)
            add_to_expression!(
                model[:obj],
                segprod[gn, t, k],
--- a/src/model/formulations/CarArr2006/structs.jl
+++ b/src/model/formulations/CarArr2006/structs.jl
@@ -14,16 +14,6 @@ module CarArr2006

 import ..PiecewiseLinearCostsFormulation

-"""
-Based on Garver (1962) and Carrión and Arryo (2006),
-which replaces (42) in Knueven et al. (2020) with a weaker version missing the on/off variable.
-Equations (45), (43), (44) in Knueven et al. (2020).
-
-Constraints
---
-* `eq_prod_above_def`: Equation (43) in Knueven et al. (2020)
-* `eq_segprod_limit`: Equation (45) in Knueven et al. (2020)
-"""
 struct PwlCosts <: PiecewiseLinearCostsFormulation end

 end
--- a/src/model/formulations/DamKucRajAta2016/ramp.jl
+++ b/src/model/formulations/DamKucRajAta2016/ramp.jl
@@ -2,26 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_ramp_eqs!
-
-Ensure constraints on ramping are met.
-Based on Damcı-Kurt et al. (2016).
-Eqns. (35), (36) in Knueven et al. (2020).
-
-Variables
---
-* :prod_above
-* :reserve
-* :is_on
-* :switch_on
-* :switch_off],
-
-Constraints
---
-* :eq_str_ramp_up
-* :eq_str_ramp_down
-"""
 function _add_ramp_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -34,14 +14,12 @@ function _add_ramp_eqs!(
    RESERVES_WHEN_RAMP_UP = true
    RESERVES_WHEN_RAMP_DOWN = true
    RESERVES_WHEN_SHUT_DOWN = true
-    is_initially_on = _is_initially_on(g)
-
-    # The following are the same for generator g across all time periods
-    SU = g.startup_limit   # startup rate
-    SD = g.shutdown_limit  # shutdown rate
-    RU = g.ramp_up_limit   # ramp up rate
-    RD = g.ramp_down_limit # ramp down rate
-
+    known_initial_conditions = true
+    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
    eq_str_ramp_down = _init(model, :eq_str_ramp_down)
    eq_str_ramp_up = _init(model, :eq_str_ramp_up)
@@ -78,7 +56,7 @@ function _add_ramp_eqs!(
        if t > 1 && time_invariant
            min_prod_last_period = prod_above[gn, t-1]

-            # Equation (35) in Knueven et al. (2020)
+            # Equation (35) in Kneuven et al. (2020)
            # Sparser version of (24)
            eq_str_ramp_up[gn, t] = @constraint(
                model,
@@ -98,7 +76,7 @@ function _add_ramp_eqs!(
            # (instead of using the amount above minimum, as min prod for t < 1 is unknown)
            max_prod_this_period += g.min_power[t] * is_on[gn, t]

-            # Modified version of equation (35) in Knueven et al. (2020)
+            # Modified version of equation (35) in Kneuven et al. (2020)
            # Equivalent to (24)
            eq_str_ramp_up[gn, t] = @constraint(
                model,
@@ -116,12 +94,12 @@ function _add_ramp_eqs!(
        on_last_period = 0.0
        if t > 1
            on_last_period = is_on[gn, t-1]
-        elseif is_initially_on
+        elseif (known_initial_conditions && g.initial_status > 0)
            on_last_period = 1.0
        end

        if t > 1 && time_invariant
-            # Equation (36) in Knueven et al. (2020)
+            # Equation (36) in Kneuven et al. (2020)
            eq_str_ramp_down[gn, t] = @constraint(
                model,
                max_prod_last_period - min_prod_this_period <=
@@ -132,7 +110,7 @@ function _add_ramp_eqs!(
            # Add back in min power
            min_prod_this_period += g.min_power[t] * is_on[gn, t]

-            # Modified version of equation (36) in Knueven et al. (2020)
+            # Modified version of equation (36) in Kneuven et al. (2020)
            # Equivalent to (25)
            eq_str_ramp_down[gn, t] = @constraint(
                model,
--- a/src/model/formulations/Gar1962/prod.jl
+++ b/src/model/formulations/Gar1962/prod.jl
@@ -2,11 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_production_vars!(model, unit, formulation_prod_vars)
-
-Adds symbols identified by `Gar1962.ProdVars` to `model`.
-"""
 function _add_production_vars!(
    model::JuMP.Model,
    g::Unit,
@@ -23,23 +18,6 @@ function _add_production_vars!(
    return
 end

-"""
-    _add_production_limit_eqs!(model, unit, formulation_prod_vars)
-
-Ensure production limit constraints are met.
-Based on Garver (1962) and Morales-España et al. (2013).
-Eqns. (18), part of (69) in Knueven et al. (2020).
-
-===
-Variables
-* :is_on
-* :prod_above
-* :reserve
-
-===
-Constraints
-* :eq_prod_limit
-"""
 function _add_production_limit_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -52,13 +30,13 @@ function _add_production_limit_eqs!(
    gn = g.name
    for t in 1:model[:instance].time
        # Objective function terms for production costs
-        # Part of (69) of Knueven et al. (2020) as C^R_g * u_g(t) term
+        # 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 Knueven et al. (2020)
+        # 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 Knueven et al. (2020)
+        # 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
--- a/src/model/formulations/Gar1962/pwlcosts.jl
+++ b/src/model/formulations/Gar1962/pwlcosts.jl
@@ -2,28 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_production_piecewise_linear_eqs!
-
-Ensure respect of production limits along each segment.
-Based on Garver (1962).
-Equations (42), (43), (44) in Knueven et al. (2020).
-NB: when reading instance, UnitCommitment.jl already calculates difference between max power for segments k and k-1,
-so the value of cost_segments[k].mw[t] is the max production *for that segment*.
-
-Added to `model`: `:eq_prod_above_def` and `:eq_segprod_limit`.
-
-===
-Variables
-* :segprod
-* :is_on
-* :prod_above
-
-===
-Constraints
-* :eq_prod_above_def
-* :eq_segprod_limit
-"""
 function _add_production_piecewise_linear_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -45,14 +23,14 @@ function _add_production_piecewise_linear_eqs!(
    K = length(g.cost_segments)
    for t in 1:model[:instance].time
        # Definition of production
-        # Equation (43) in Knueven et al. (2020)
+        # Equation (43) in Kneuven et al. (2020)
        eq_prod_above_def[gn, t] = @constraint(
            model,
            prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
        )

        for k in 1:K
-            # Equation (42) in Knueven et al. (2020)
+            # Equation (42) in Kneuven et al. (2020)
            # Without this, solvers will add a lot of implied bound cuts to
            # have this same effect.
            # NB: when reading instance, UnitCommitment.jl already calculates
@@ -69,7 +47,7 @@ function _add_production_piecewise_linear_eqs!(
            set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])

            # Objective function
-            # Equation (44) in Knueven et al. (2020)
+            # Equation (44) in Kneuven et al. (2020)
            add_to_expression!(
                model[:obj],
                segprod[gn, t, k],
--- a/src/model/formulations/Gar1962/status.jl
+++ b/src/model/formulations/Gar1962/status.jl
@@ -2,12 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_status_vars!
-
-Adds symbols identified by `Gar1962.StatusVars` to `model`.
-Fix variables if a certain generator _must_ run or based on initial conditions.
-"""
 function _add_status_vars!(
    model::JuMP.Model,
    g::Unit,
@@ -16,108 +10,20 @@ function _add_status_vars!(
    is_on = _init(model, :is_on)
    switch_on = _init(model, :switch_on)
    switch_off = _init(model, :switch_off)
-    FIX_VARS = !formulation_status_vars.fix_vars_via_constraint
-    is_initially_on = _is_initially_on(g) > 0
    for t in 1:model[:instance].time
-        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)
-
-        # Use initial conditions and whether a unit must run to fix variables
-        if FIX_VARS
-            # Fix variables using fix function
-            if g.must_run[t]
-                # If the generator _must_ run, then it is obviously on and cannot be switched off
-                # In the first time period, force unit to switch on if was off before
-                # Otherwise, unit is on, and will never turn off, so will never need to turn on
-                fix(is_on[g.name, t], 1.0; force = true)
-                fix(
-                    switch_on[g.name, t],
-                    (t == 1 ? 1.0 - _is_initially_on(g) : 0.0);
-                    force = true,
-                )
-                fix(switch_off[g.name, t], 0.0; force = true)
-            elseif t == 1
-                if is_initially_on
-                    # Generator was on (for g.initial_status time periods),
-                    # so cannot be more switched on until the period after the first time it can be turned off
-                    fix(switch_on[g.name, 1], 0.0; force = true)
-                else
-                    # Generator is initially off (for -g.initial_status time periods)
-                    # Cannot be switched off more
-                    fix(switch_off[g.name, 1], 0.0; force = true)
-                end
-            end
+        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
-            # Add explicit constraint if !FIX_VARS
-            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
-            elseif t == 1
-                if is_initially_on
-                    switch_on[g.name, t] = 0.0
-                else
-                    switch_off[g.name, t] = 0.0
-                end
-            end
-        end
-
-        # Use initial conditions and whether a unit must run to fix variables
-        if FIX_VARS
-            # Fix variables using fix function
-            if g.must_run[t]
-                # If the generator _must_ run, then it is obviously on and cannot be switched off
-                # In the first time period, force unit to switch on if was off before
-                # Otherwise, unit is on, and will never turn off, so will never need to turn on
-                fix(is_on[g.name, t], 1.0; force = true)
-                fix(
-                    switch_on[g.name, t],
-                    (t == 1 ? 1.0 - _is_initially_on(g) : 0.0);
-                    force = true,
-                )
-                fix(switch_off[g.name, t], 0.0; force = true)
-            elseif t == 1
-                if is_initially_on
-                    # Generator was on (for g.initial_status time periods),
-                    # so cannot be more switched on until the period after the first time it can be turned off
-                    fix(switch_on[g.name, 1], 0.0; force = true)
-                else
-                    # Generator is initially off (for -g.initial_status time periods)
-                    # Cannot be switched off more
-                    fix(switch_off[g.name, 1], 0.0; force = true)
-                end
-            end
-        else
-            # Add explicit constraint if !FIX_VARS
-            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
-            elseif t == 1
-                if is_initially_on
-                    switch_on[g.name, t] = 0.0
-                else
-                    switch_off[g.name, t] = 0.0
-                end
-            end
+            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

-"""
-    _add_status_eqs!
-
-Creates constraints `eq_binary_link` and `eq_switch_on_off` using variables in `Gar1962.StatusVars`.
-
-Constraints
---
-* `eq_binary_link`
-* `eq_switch_on_off`
-"""
 function _add_status_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -129,32 +35,27 @@ function _add_status_eqs!(
    switch_off = model[:switch_off]
    switch_on = model[:switch_on]
    for t in 1:model[:instance].time
-        if g.must_run[t]
-            continue
-        end
-
-        # Link binary variables
-        # Equation (2) in Knueven et al. (2020), originally from Garver (1962)
-        if t == 1
-            eq_binary_link[g.name, t] = @constraint(
+        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,
-                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]
+                switch_on[g.name, t] + switch_off[g.name, t] <= 1
            )
        end
-
-        # Cannot switch on and off at the same time
-        # amk: I am not sure this is in Knueven et al. (2020)
-        eq_switch_on_off[g.name, t] = @constraint(
-            model,
-            switch_on[g.name, t] + switch_off[g.name, t] <= 1
-        )
    end
    return
 end
--- a/src/model/formulations/Gar1962/structs.jl
+++ b/src/model/formulations/Gar1962/structs.jl
@@ -17,53 +17,8 @@ import ..PiecewiseLinearCostsFormulation
 import ..ProductionVarsFormulation
 import ..StatusVarsFormulation

-"""
-Variables
---
-* `prod_above`:
-        [gen, t];
-        *production above minimum required level*;
-        lb: 0, ub: Inf.
-        KnuOstWat2020: `p'_g(t)`
-* `segprod`:
-        [gen, segment, t];
-        *how much generator produces on cost segment in time t*;
-        lb: 0, ub: Inf.
-        KnuOstWat2020: `p_g^l(t)`
-"""
 struct ProdVars <: ProductionVarsFormulation end
-
 struct PwlCosts <: PiecewiseLinearCostsFormulation end
-
-"""
-Variables
---
-* `is_on`:
-        [gen, t];
-        *is generator on at time t?*
-        lb: 0, ub: 1, binary.
-        KnuOstWat2020: `u_g(t)`
-* `switch_on`:
-        [gen, t];
-        *indicator that generator will be turned on at t*;
-        lb: 0, ub: 1, binary.
-        KnuOstWat2020: `v_g(t)`
-* `switch_off`: binary;
-        [gen, t];
-        *indicator that generator will be turned off at t*;
-        lb: 0, ub: 1, binary.
-        KnuOstWat2020: `w_g(t)`
-
-Arguments
---
-* `fix_vars_via_constraint`:
-        indicator for whether to set vars to a constant using `fix` or by adding an explicit constraint
-        (particulary useful for debugging purposes).
-"""
-struct StatusVars <: StatusVarsFormulation
-    fix_vars_via_constraint::Bool
-
-    StatusVars() = new(false)
-end
+struct StatusVars <: StatusVarsFormulation end

 end
--- a/src/model/formulations/GenMorRam2017/startstop.jl
+++ b/src/model/formulations/GenMorRam2017/startstop.jl
@@ -1,100 +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.
-
-"""
-    _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
-
-Startup and shutdown limits from Gentile et al. (2017).
-Eqns. (20), (23a), and (23b) in Knueven et al. (2020).
-
-Creates constraints `eq_startstop_limit`,  `eq_startup_limit`, and `eq_shutdown_limit`
-using variables `Gar1962.StatusVars`, `prod_above` from `Gar1962.ProdVars`, and `reserve`.
-
-Constraints
---
-* `eq_startstop_limit`
-* `eq_startup_limit`
-* `eq_shutdown_limit`
-"""
-function _add_startup_shutdown_limit_eqs!(
-    model::JuMP.Model,
-    g::Unit,
-    formulation_status_vars::Gar1962.StatusVars,
-)::Nothing
-    # TODO: Move upper case constants to model[:instance]
-    RESERVES_WHEN_START_UP = true
-    RESERVES_WHEN_RAMP_UP = true
-    RESERVES_WHEN_RAMP_DOWN = true
-    RESERVES_WHEN_SHUT_DOWN = true
-
-    eq_startstop_limit = _init(model, :eq_startstop_limit)
-    eq_shutdown_limit = _init(model, :eq_shutdown_limit)
-    eq_startup_limit = _init(model, :eq_startup_limit)
-
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
-    reserve = model[:reserve]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
-
-    T = model[:instance].time
-    gi = g.name
-
-    if g.initial_power > g.shutdown_limit
-        #eqs.shutdown_limit[gi, 0] = @constraint(mip, vars.switch_off[gi, 1] <= 0)
-        if formulation_status_vars.always_create_vars
-            fix(switch_off[gi, 1], 0.0; force = true)
-            @constraint(mip, vars.switch_off[gi, 1] <= 0)
-        else
-            switch_off[gi, 1] = 0.0
-        end
-    end
-
-    for t in 1:T
-        ## 2020-10-09 amk: added eqn (20) and check of g.min_uptime
-        # Not present in (23) in Kneueven et al.
-        if g.min_uptime > 1
-            # Equation (20) in Knueven et al. (2020)
-            eqs.startstop_limit[gi, t] = @constraint(
-                model,
-                prod_above[gi, t] + reserve[gi, t] <=
-                (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
-                max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t] - (
-                    t < T ?
-                    max(0, g.max_power[t] - g.shutdown_limit) *
-                    switch_off[gi, t+1] : 0.0
-                )
-            )
-        else
-            ## Startup limits
-            # Equation (23a) in Knueven et al. (2020)
-            eqs.startup_limit[gi, t] = @constraint(
-                model,
-                prod_above[gi, t] + reserve[gi, t] <=
-                (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
-                max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t] - (
-                    t < T ?
-                    max(0, g.startup_limit - g.shutdown_limit) *
-                    switch_off[gi, t+1] : 0.0
-                )
-            )
-
-            ## Shutdown limits
-            if t < T
-                # Equation (23b) in Knueven et al. (2020)
-                eqs.shutdown_limit[gi, t] = @constraint(
-                    model,
-                    prod_above[gi, t] + reserve[gi, t] <=
-                    (g.max_power[t] - g.min_power[t]) * xis_on[gi, t] - (
-                        t < T ?
-                        max(0, g.max_power[t] - g.shutdown_limit) *
-                        switch_off[gi, t+1] : 0.0
-                    ) -
-                    max(0, g.shutdown_limit - g.startup_limit) *
-                    switch_on[gi, t]
-                )
-            end
-        end
-    end
-end
--- a/src/model/formulations/KnuOstWat2018/pwlcosts.jl
+++ b/src/model/formulations/KnuOstWat2018/pwlcosts.jl
@@ -2,30 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_production_piecewise_linear_eqs!
-
-Ensure respect of production limits along each segment.
-Based on Knueven et al. (2018b).
-Eqns. (43), (44), (46), (48) in Knueven et al. (2020).
-NB: when reading instance, UnitCommitment.jl already calculates difference between max power for segments k and k-1
-so the value of cost_segments[k].mw[t] is the max production *for that segment*.
-
-===
-Variables
-* :segprod
-* :is_on
-* :switch_on
-* :switch_off
-* :prod_above
-
-===
-Constraints
-* :eq_prod_above_def
-* :eq_segprod_limit_a
-* :eq_segprod_limit_b
-* :eq_segprod_limit_c
-"""
 function _add_production_piecewise_linear_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -81,7 +57,7 @@ function _add_production_piecewise_linear_eqs!(
            end

            if g.min_uptime > 1
-                # Equation (46) in Knueven et al. (2020)
+                # Equation (46) in Kneuven et al. (2020)
                eq_segprod_limit_a[gn, t, k] = @constraint(
                    model,
                    segprod[gn, t, k] <=
@@ -90,7 +66,7 @@ function _add_production_piecewise_linear_eqs!(
                    (t < T ? Cw * switch_off[gn, t+1] : 0.0)
                )
            else
-                # Equation (47a)/(48a) in Knueven et al. (2020)
+                # Equation (47a)/(48a) in Kneuven et al. (2020)
                eq_segprod_limit_b[gn, t, k] = @constraint(
                    model,
                    segprod[gn, t, k] <=
@@ -99,7 +75,7 @@ function _add_production_piecewise_linear_eqs!(
                    (t < T ? max(0, Cv - Cw) * switch_off[gn, t+1] : 0.0)
                )

-                # Equation (47b)/(48b) in Knueven et al. (2020)
+                # Equation (47b)/(48b) in Kneuven et al. (2020)
                eq_segprod_limit_c[gn, t, k] = @constraint(
                    model,
                    segprod[gn, t, k] <=
@@ -110,14 +86,14 @@ function _add_production_piecewise_linear_eqs!(
            end

            # Definition of production
-            # Equation (43) in Knueven et al. (2020)
+            # Equation (43) in Kneuven et al. (2020)
            eq_prod_above_def[gn, t] = @constraint(
                model,
                prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
            )

            # Objective function
-            # Equation (44) in Knueven et al. (2020)
+            # Equation (44) in Kneuven et al. (2020)
            add_to_expression!(
                model[:obj],
                segprod[gn, t, k],
--- a/src/model/formulations/KnuOstWat2018/scosts.jl
+++ b/src/model/formulations/KnuOstWat2018/scosts.jl
@@ -1,123 +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.
-
-"""
-    _add_startup_cost_eqs!
-
-Extended formulation of startup costs using indicator variables
-based on Knueven, Ostrowski, and Watson, 2020
--- equations (59), (60), (61).
-
-Variables
---
-* switch_on
-* switch_off
-* downtime_arc
-
-Constraints
---
-* eq_startup_at_t
-* eq_shutdown_at_t
-"""
-function _add_startup_cost_eqs!(
-    model::JuMP.Model,
-    g::Unit,
-    formulation::MorLatRam2013.StartupCosts,
-)::Nothing
-    S = length(g.startup_categories)
-    if S == 0
-        return
-    end
-    gn = g.name
-
-    _init(model, eq_startup_at_t)
-    _init(model, eq_shutdown_at_t)
-
-    switch_on = model[:switch_on]
-    switch_off = model[:switch_off]
-    downtime_arc = model[:downtime_arc]
-
-    DT = g.min_downtime # minimum time offline
-    TC = g.startup_categories[S].delay # time offline until totally cold
-
-    # If initial_status < 0, then this is the amount of time the generator has been off
-    initial_time_shutdown = (g.initial_status < 0 ? -g.initial_status : 0)
-
-    for t in 1:model[:instance].time
-        # Fix to zero values of downtime_arc outside the feasible time pairs
-        # Specifically, x(t,t') = 0 if t' does not belong to 𝒢 = [t+DT, t+TC-1]
-        # This is because DT is the minimum downtime, so there is no way x(t,t')=1 for t'<t+DT
-        # and TC is the "time until cold" => if the generator starts afterwards, always has max cost
-        #start_time = min(t + DT, T)
-        #end_time = min(t + TC - 1, T)
-        #for tmp_t in t+1:start_time
-        #  fix(vars.downtime_arc[gn, t, tmp_t], 0.; force = true)
-        #end
-        #for tmp_t in end_time+1:T
-        #  fix(vars.downtime_arc[gn, t, tmp_t], 0.; force = true)
-        #end
-
-        # Equation (59) in Knueven et al. (2020)
-        # Relate downtime_arc with switch_on
-        # "switch_on[g,t] >= x_g(t',t) for all t' \in [t-TC+1, t-DT]"
-        eq_startup_at_t[gn, t] = @constraint(
-            model,
-            switch_on[gn, t] >= sum(
-                downtime_arc[gn, tmp_t, t] for
-                tmp_t in t-TC+1:t-DT if tmp_t >= 1
-            )
-        )
-
-        # Equation (60) in Knueven et al. (2020)
-        # "switch_off[g,t] >= x_g(t,t') for all t' \in [t+DT, t+TC-1]"
-        eqs.shutdown_at_t[gn, t] = @constraint(
-            model,
-            switch_off[gn, t] >= sum(
-                downtime_arc[gn, t, tmp_t] for
-                tmp_t in t+DT:t+TC-1 if tmp_t <= T
-            )
-        )
-
-        # Objective function terms for start-up costs
-        # Equation (61) in Knueven et al. (2020)
-        default_category = S
-        if initial_time_shutdown > 0 && t + initial_time_shutdown - 1 < TC
-            for s in 1:S-1
-                # If off for x periods before, then belongs to category s
-                # if -x+1 in [t-delay[s+1]+1,t-delay[s]]
-                # or, equivalently, if total time off in [delay[s], delay[s+1]-1]
-                # where total time off = t - 1 + initial_time_shutdown
-                # (the -1 because not off for current time period)
-                if t + initial_time_shutdown - 1 <
-                   g.startup_categories[s+1].delay
-                    default_category = s
-                    break # does not go into next category
-                end
-            end
-        end
-        add_to_expression!(
-            model[:obj],
-            switch_on[gn, t],
-            g.startup_categories[default_category].cost,
-        )
-
-        for s in 1:S-1
-            # Objective function terms for start-up costs
-            # Equation (61) in Knueven et al. (2020)
-            # Says to replace the cost of last category with cost of category s
-            start_range = max((t - g.startup_categories[s+1].delay + 1), 1)
-            end_range = min((t - g.startup_categories[s].delay), T - 1)
-            for tmp_t in start_range:end_range
-                if (t < tmp_t + DT) || (t >= tmp_t + TC) # the second clause should never be true for s < S
-                    continue
-                end
-                add_to_expression!(
-                    model[:obj],
-                    downtime_arc[gn, tmp_t, t],
-                    g.startup_categories[s].cost - g.startup_categories[S].cost,
-                )
-            end
-        end # iterate over startup categories
-    end # iterate over time
-end # add_startup_costs_KneOstWat20
--- a/src/model/formulations/MorLatRam2013/ramp.jl
+++ b/src/model/formulations/MorLatRam2013/ramp.jl
@@ -2,27 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_ramp_eqs!
-
-Ensure constraints on ramping are met.
-Needs to be used in combination with shutdown rate constraints, e.g., (21b) in Knueven et al. (2020).
-Based on Morales-España, Latorre, and Ramos, 2013.
-Eqns. (26)+(27) [replaced by (24)+(25) if time-varying min demand] in Knueven et al. (2020).
-
-Variables
---
-* :is_on
-* :switch_off
-* :switch_on
-* :prod_above
-* :reserve
-
-Constraints
---
-* :eq_ramp_up
-* :eq_ramp_down
-"""
 function _add_ramp_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -36,10 +15,10 @@ function _add_ramp_eqs!(
    RESERVES_WHEN_RAMP_DOWN = true
    RESERVES_WHEN_SHUT_DOWN = true
    is_initially_on = (g.initial_status > 0)
-    SU = g.startup_limit   # startup rate
-    SD = g.shutdown_limit  # shutdown rate
-    RU = g.ramp_up_limit   # ramp up rate
-    RD = g.ramp_down_limit # ramp down rate
+    SU = g.startup_limit
+    SD = g.shutdown_limit
+    RU = g.ramp_up_limit
+    RD = g.ramp_down_limit
    gn = g.name
    eq_ramp_down = _init(model, :eq_ramp_down)
    eq_ramp_up = _init(model, :eq_str_ramp_up)
@@ -92,7 +71,7 @@ function _add_ramp_eqs!(
                    RU * is_on[gn, t-1] + SU * switch_on[gn, t]
                )
            else
-                # Equation (26) in Knueven et al. (2020)
+                # Equation (26) in Kneuven et al. (2020)
                # TODO: what if RU < SU? places too stringent upper bound
                # prod_above[gn, t] when starting up, and creates diff with (24).
                eq_ramp_up[gn, t] = @constraint(
@@ -136,7 +115,7 @@ function _add_ramp_eqs!(
                    RD * is_on[gn, t] + SD * switch_off[gn, t]
                )
            else
-                # Equation (27) in Knueven et al. (2020)
+                # Equation (27) in Kneuven et al. (2020)
                # TODO: Similar to above, what to do if shutting down in time t
                # and RD < SD? There is a difference with (25).
                eq_ramp_down[gn, t] = @constraint(
--- a/src/model/formulations/MorLatRam2013/scosts.jl
+++ b/src/model/formulations/MorLatRam2013/scosts.jl
@@ -2,56 +2,21 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_startup_cost_eqs!
-
-Extended formulation of startup costs using indicator variables
-based on Muckstadt and Wilson, 1968;
-this version by Morales-España, Latorre, and Ramos, 2013.
-Eqns. (54), (55), and (56) in Knueven et al. (2020).
-Note that the last 'constraint' is actually setting the objective.
-
-\tstartup[gi,s,t] ≤ sum_{i=s.delay}^{(s+1).delay-1} switch_off[gi,t-i]
-\tswitch_on[gi,t] = sum_{s=1}^{length(startup_categories)} startup[gi,s,t]
-\tstartup_cost[gi,t] = sum_{s=1}^{length(startup_categories)} cost_segments[s].cost * startup[gi,s,t]
-
-Variables
---
-* startup
-* switch_on
-* switch_off
-
-Constraints
---
-* eq_startup_choose
-* eq_startup_restrict
-"""
 function _add_startup_cost_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation::MorLatRam2013.StartupCosts,
 )::Nothing
-    S = length(g.startup_categories)
-    if S == 0
-        return
-    end
-
-    # Constraints created
    eq_startup_choose = _init(model, :eq_startup_choose)
    eq_startup_restrict = _init(model, :eq_startup_restrict)
-
-    # Variables needed
+    S = length(g.startup_categories)
    startup = model[:startup]
-    switch_on = model[:switch_on]
-    switch_off = model[:switch_off]
-
-    gn = g.name
    for t in 1:model[:instance].time
        # If unit is switching on, we must choose a startup category
-        # Equation (55) in Knueven et al. (2020)
-        eq_startup_choose[gn, t] = @constraint(
+        eq_startup_choose[g.name, t] = @constraint(
            model,
-            switch_on[gn, t] == sum(startup[gn, t, s] for s in 1:S)
+            model[:switch_on][g.name, t] ==
+            sum(startup[g.name, t, s] for s in 1:S)
        )

        for s in 1:S
@@ -61,28 +26,25 @@ function _add_startup_cost_eqs!(
                range_start = t - g.startup_categories[s+1].delay + 1
                range_end = t - g.startup_categories[s].delay
                range = (range_start:range_end)
-                # If initial_status < 0, then this is the amount of time the generator has been off
                initial_sum = (
                    g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0
                )
-                # Change of index version of equation (54) in Knueven et al. (2020):
-                #   startup[gi,s,t] ≤ sum_{i=s.delay}^{(s+1).delay-1} switch_off[gi,t-i]
-                eq_startup_restrict[gn, t, s] = @constraint(
+                eq_startup_restrict[g.name, t, s] = @constraint(
                    model,
-                    startup[gn, t, s] <=
-                    initial_sum +
-                    sum(switch_off[gn, i] for i in range if i >= 1)
+                    startup[g.name, t, s] <=
+                    initial_sum + sum(
+                        model[:switch_off][g.name, i] for i in range if i >= 1
+                    )
                )
-            end # if s < S (not the last category)
+            end

            # Objective function terms for start-up costs
-            # Equation (56) in Knueven et al. (2020)
            add_to_expression!(
                model[:obj],
-                startup[gn, t, s],
+                startup[g.name, t, s],
                g.startup_categories[s].cost,
            )
-        end # iterate over startup categories
-    end # iterate over time
+        end
+    end
    return
 end
--- a/src/model/formulations/MorLatRam2013/startstop.jl
+++ b/src/model/formulations/MorLatRam2013/startstop.jl
@@ -1,89 +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.
-
-"""
-    _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
-
-Startup and shutdown limits from Morales-España et al. (2013a).
-Eqns. (20), (21a), and (21b) in Knueven et al. (2020).
-
-Variables
---
-* :is_on
-* :prod_above
-* :reserve
-* :switch_on
-* :switch_off
-
-Constraints
---
-* :eq_startstop_limit
-* :eq_startup_limit
-* :eq_shutdown_limit
-"""
-function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    # TODO: Move upper case constants to model[:instance]
-    RESERVES_WHEN_START_UP = true
-    RESERVES_WHEN_RAMP_UP = true
-    RESERVES_WHEN_RAMP_DOWN = true
-    RESERVES_WHEN_SHUT_DOWN = true
-
-    eq_startstop_limit = _init(model, :eq_startstop_limit)
-    eq_shutdown_limit = _init(model, :eq_shutdown_limit)
-    eq_startup_limit = _init(model, :eq_startup_limit)
-
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
-    reserve = model[:reserve]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
-
-    T = model[:instance].time
-    gi = g.name
-    for t in 1:T
-        ## 2020-10-09 amk: added eqn (20) and check of g.min_uptime
-        if g.min_uptime > 1 && t < T
-            # Equation (20) in Knueven et al. (2020)
-            # UT > 1 required, to guarantee that vars.switch_on[gi, t] and vars.switch_off[gi, t+1] are not both = 1 at the same time
-            eq_startstop_limit[gi, t] = @constraint(
-                model,
-                prod_above[gi, t] + reserve[gi, t] <=
-                (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
-                max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t] -
-                max(0, g.max_power[t] - g.shutdown_limit) * switch_off[gi, t+1]
-            )
-        else
-            ## Startup limits
-            # Equation (21a) in Knueven et al. (2020)
-            # Proposed by Morales-España et al. (2013a)
-            eqs_startup_limit[gi, t] = @constraint(
-                model,
-                prod_above[gi, t] + reserve[gi, t] <=
-                (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
-                max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t]
-            )
-
-            ## Shutdown limits
-            if t < T
-                # Equation (21b) in Knueven et al. (2020)
-                # TODO different from what was in previous model, due to reserve variable
-                # ax: ideally should have reserve_up and reserve_down variables
-                #     i.e., the generator should be able to increase/decrease production as specified
-                #     (this is a heuristic for a "robust" solution,
-                #     in case there is an outage or a surge, and flow has to be redirected)
-                # amk: if shutdown_limit is the max prod of generator in time period before shutting down,
-                #      then it makes sense to count reserves, because otherwise, if reserves ≠ 0,
-                #      then the generator will actually produce more than the limit
-                eqs.shutdown_limit[gi, t] = @constraint(
-                    model,
-                    prod_above[gi, t] +
-                    (RESERVES_WHEN_SHUT_DOWN ? reserve[gi, t] : 0.0) <=
-                    (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
-                    max(0, g.max_power[t] - g.shutdown_limit) *
-                    switch_off[gi, t+1]
-                )
-            end
-        end # check if g.min_uptime > 1
-    end # loop over time
-end # _add_startup_shutdown_limit_eqs!
--- a/src/model/formulations/NowRom2000/scosts.jl
+++ b/src/model/formulations/NowRom2000/scosts.jl
@@ -1,20 +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.
-
-"""
-    _add_startup_cost_eqs!
-
-Based on Nowak and Römisch, 2000.
-Introduces auxiliary startup cost variable, c_g^SU(t) for each time period,
-and uses startup status variable, u_g(t);
-there are exponentially many facets in this space,
-but there is a linear-time separation algorithm (Brandenburg et al., 2017).
-"""
-function _add_startup_cost_eqs!(
-    model::JuMP.Model,
-    g::Unit,
-    formulation::MorLatRam2013.StartupCosts,
-)::Nothing
-    return error("Not implemented.")
-end
--- a/src/model/formulations/OstAnjVan2012/ramp.jl
+++ b/src/model/formulations/OstAnjVan2012/ramp.jl
@@ -1,85 +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.
-
-"""
-    _add_ramp_eqs!
-
-Ensure constraints on ramping are met.
-Based on Ostrowski, Anjos, Vannelli (2012).
-Eqn (37) in Knueven et al. (2020).
-
-Variables
---
-* :is_on
-* :prod_above
-* :reserve
-
-Constraints
---
-* :eq_str_prod_limit
-"""
-function _add_ramp_eqs!(
-    model::JuMP.Model,
-    g::Unit,
-    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
-    RESERVES_WHEN_RAMP_UP = true
-    RESERVES_WHEN_RAMP_DOWN = true
-    RESERVES_WHEN_SHUT_DOWN = true
-    is_initially_on = _is_initially_on(g)
-
-    gn = g.name
-    eq_str_prod_limit = _init(model, :eq_str_prod_limit)
-
-    # Variables that we need
-    reserve = model[:reserve]
-
-    # Gar1962.ProdVars
-    prod_above = model[:prod_above]
-
-    # Gar1962.StatusVars
-    is_on = model[:is_on]
-    switch_off = model[:switch_off]
-
-    # The following are the same for generator g across all time periods
-    UT = g.min_uptime
-
-    SU = g.startup_limit   # startup rate
-    SD = g.shutdown_limit  # shutdown rate
-    RU = g.ramp_up_limit   # ramp up rate
-    RD = g.ramp_down_limit # ramp down rate
-
-    # TODO check initial conditions, but maybe okay as long as (35) and (36) are also used
-    for t in 1:model[:instance].time
-        Pbar = g.max_power[t]
-
-        #TRD = floor((Pbar - SU)/RD)
-        # TODO check amk changed TRD wrt Knueven et al.
-        TRD = ceil((Pbar - SD) / RD) # ramp down time
-
-        if Pbar < 1e-7
-            # Skip this time period if max power = 0
-            continue
-        end
-
-        if UT >= 1
-            # Equation (37) in Knueven et al. (2020)
-            KSD = min(TRD, UT - 1, T - t - 1)
-            eq_str_prod_limit[gn, t] = @constraint(
-                model,
-                prod_above[gn, t] +
-                g.min_power[t] * is_on[gn, t] +
-                (RESERVES_WHEN_RAMP_DOWN ? reserve[gn, t] : 0.0) <=
-                Pbar * is_on[gi, t] - sum(
-                    (Pbar - (SD + i * RD)) * switch_off[gi, t+1+i] for
-                    i in 0:KSD
-                )
-            )
-        end # check UT >= 1
-    end # loop over time
-end
--- a/src/model/formulations/PanGua2016/ramp.jl
+++ b/src/model/formulations/PanGua2016/ramp.jl
@@ -2,28 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_ramp_eqs!
-
-Add tighter upper bounds on production based on ramp-down trajectory.
-Based on (28) in Pan and Guan (2016).
-But there is an extra time period covered using (40) of Knueven et al. (2020).
-Eqns. (38), (40), (41) in Knueven et al. (2020).
-
-Variables
---
-* :prod_above
-* :reserve
-* :is_on
-* :switch_on
-* :switch_off
-
-Constraints
---
-* :str_prod_limit
-* :prod_limit_ramp_up_extra_period
-* :prod_limit_shutdown_trajectory
-"""
 function _add_ramp_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -63,14 +41,14 @@ function _add_ramp_eqs!(
        end

        #TRD = floor((Pbar - SU) / RD) # ramp down time
-        # TODO check amk changed TRD wrt Knueven et al.
+        # TODO check amk changed TRD wrt Kneuven et al.
        TRD = ceil((Pbar - SD) / RD)  # ramp down time
        TRU = floor((Pbar - SU) / RU) # ramp up time, can be negative if Pbar < SU

        # TODO check initial time periods: what if generator has been running for x periods?
        # But maybe ok as long as (35) and (36) are also used...
        if UT > 1
-            # Equation (38) in Knueven et al. (2020)
+            # Equation (38) in Kneuven et al. (2020)
            # Generalization of (20)
            # Necessary that if any of the switch_on = 1 in the sum,
            # then switch_off[gn, t+1] = 0
@@ -87,7 +65,7 @@ function _add_ramp_eqs!(
            )

            if UT - 2 < TRU
-                # Equation (40) in Knueven et al. (2020)
+                # Equation (40) in Kneuven et al. (2020)
                # Covers an additional time period of the ramp-up trajectory, compared to (38)
                eq_prod_limit_ramp_up_extra_period[gn, t] = @constraint(
                    model,
@@ -105,7 +83,7 @@ function _add_ramp_eqs!(
            KSD = min(TRD, UT - 1, T - t - 1)
            if KSD > 0
                KSU = min(TRU, UT - 2 - KSD, t - 1)
-                # Equation (41) in Knueven et al. (2020)
+                # Equation (41) in Kneuven et al. (2020)
                eq_prod_limit_shutdown_trajectory[gn, t] = @constraint(
                    model,
                    prod_above[gn, t] +
--- a/src/model/formulations/base/bus.jl
+++ b/src/model/formulations/base/bus.jl
@@ -2,11 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_bus!(model::JuMP.Model, b::Bus)::Nothing
-
-Creates `expr_net_injection` and adds `curtail` variable to `model`.
-"""
 function _add_bus!(model::JuMP.Model, b::Bus)::Nothing
    net_injection = _init(model, :expr_net_injection)
    curtail = _init(model, :curtail)
--- a/src/model/formulations/base/psload.jl
+++ b/src/model/formulations/base/psload.jl
@@ -2,9 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_price_sensitive_load!(model::JuMP.Model, ps::PriceSensitiveLoad)
-"""
 function _add_price_sensitive_load!(
    model::JuMP.Model,
    ps::PriceSensitiveLoad,
--- a/src/model/formulations/base/sensitivity.jl
+++ b/src/model/formulations/base/sensitivity.jl
@@ -18,7 +18,7 @@ function _injection_shift_factors(;
    lines::Array{TransmissionLine},
 )
    susceptance = _susceptance_matrix(lines)
-    incidence = _reduced_incidence_matrix(buses = buses, lines = lines)
+    incidence = _reduced_incidence_matrix(lines = lines, buses = buses)
    laplacian = transpose(incidence) * susceptance * incidence
    isf = susceptance * incidence * inv(Array(laplacian))
    return isf
--- a/src/model/formulations/base/structs.jl
+++ b/src/model/formulations/base/structs.jl
@@ -9,15 +9,6 @@ abstract type StartupCostsFormulation end
 abstract type StatusVarsFormulation end
 abstract type ProductionVarsFormulation end

-"""
-    struct Formulation
-
-Every formulation has to specify components for setting production variables and limits,
-startup costs and piecewise-linear costs, ramping variables and constraints,
-status variables (on/off, starting up, shutting down), and transmission constraints.
-
-Some of these components are allowed to be empty, as long as overall validity of the formulation is maintained.
-"""
 struct Formulation
    prod_vars::ProductionVarsFormulation
    pwl_costs::PiecewiseLinearCostsFormulation
--- a/src/model/formulations/base/system.jl
+++ b/src/model/formulations/base/system.jl
@@ -2,32 +2,12 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_system_wide_eqs!(model::JuMP.Model)::Nothing
-
-Adds constraints that apply to the whole system, such as relating to net injection and reserves.
-"""
 function _add_system_wide_eqs!(model::JuMP.Model)::Nothing
    _add_net_injection_eqs!(model)
    _add_reserve_eqs!(model)
    return
 end

-"""
-    _add_net_injection_eqs!(model::JuMP.Model)::Nothing
-
-Adds `net_injection`, `eq_net_injection_def`, and `eq_power_balance` identifiers into `model`.
-
-Variables
---
-* `expr_net_injection`
-* `net_injection`
-
-Constraints
---
-* `eq_net_injection_def`
-* `eq_power_balance`
-"""
 function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
    T = model[:instance].time
    net_injection = _init(model, :net_injection)
@@ -47,30 +27,11 @@ function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
    return
 end

-"""
-    _add_reserve_eqs!(model::JuMP.Model)::Nothing
-
-Ensure constraints on reserves are met.
-Based on Morales-España et al. (2013a).
-Eqn. (68) from Knueven et al. (2020).
-
-Adds `eq_min_reserve` identifier to `model`, and corresponding constraint.
-
-Variables
---
-* `reserve`
-* `reserve_shortfall`
-
-Constraints
---
-* `eq_min_reserve`
-"""
 function _add_reserve_eqs!(model::JuMP.Model)::Nothing
-    instance = model[:instance]
    eq_min_reserve = _init(model, :eq_min_reserve)
    instance = model[:instance]
    for t in 1:instance.time
-        # Equation (68) in Knueven et al. (2020)
+        # 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)
--- a/src/model/formulations/base/unit.jl
+++ b/src/model/formulations/base/unit.jl
@@ -2,15 +2,6 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.

-"""
-    _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
-
-Add production, reserve, startup, shutdown, and status variables,
-and constraints for min uptime/downtime, net injection, production, ramping, startup, shutdown, and status.
-
-Fix variables if a certain generator _must_ run or if a generator provides spinning reserves.
-Also, add overflow penalty to objective for each transmission line.
-"""
 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")
@@ -44,41 +35,21 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
        formulation.status_vars,
    )
    _add_startup_cost_eqs!(model, g, formulation.startup_costs)
-    _add_shutdown_cost_eqs!(model, g)
-    _add_startup_shutdown_limit_eqs!(
-        model,
-        g,
-        formulation.status_vars,
-        formulation.prod_vars,
-    )
+    _add_startup_shutdown_limit_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)

-"""
-    _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
-
-Add `:reserve` variable to `model`, fixed to zero if no spinning reserves specified.
-"""
-function _add_reserve_vars!(
-    model::JuMP.Model,
-    g::Unit,
-    ALWAYS_CREATE_VARS = false,
-)::Nothing
+function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
    reserve = _init(model, :reserve)
-    reserve_shortfall = _init(model, :reserve_shortfall) # for accounting for shortfall penalty in the objective
+    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
-            if ALWAYS_CREATE_VARS
-                reserve[g.name, t] = @variable(model, lower_bound = 0)
-                fix(reserve[g.name, t], 0.0; force = true)
-            else
-                reserve[g.name, t] = 0.0
-            end
+            reserve[g.name, t] = 0.0
        end
        reserve_shortfall[t] =
            (model[:instance].shortfall_penalty[t] >= 0) ?
@@ -87,19 +58,14 @@ function _add_reserve_vars!(
    return
 end

-"""
-    _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
-"""
 function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    # nothing to do here
+    reserve = model[:reserve]
+    for t in 1:model[:instance].time
+        add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0)
+    end
    return
 end

-"""
-    _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
-
-Add `startup` to model.
-"""
 function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
    startup = _init(model, :startup)
    for t in 1:model[:instance].time
@@ -110,31 +76,14 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
    return
 end

-"""
-    _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
-
-Creates startup/shutdown limit constraints below based on variables `Gar1962.StatusVars`, `prod_above` from `Gar1962.ProdVars`, and `reserve`.
-
-Constraints
---
-* :eq_startup_limit
-* :eq_shutdown_limit
-"""
-function _add_startup_shutdown_limit_eqs!(
-    model::JuMP.Model,
-    g::Unit,
-    formulation_status_vars::Gar1962.StatusVars,
-    formulation_prod_vars::Gar1962.ProdVars,
-)::Nothing
+function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
    eq_shutdown_limit = _init(model, :eq_shutdown_limit)
    eq_startup_limit = _init(model, :eq_startup_limit)
-
    is_on = model[:is_on]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
    prod_above = model[:prod_above]
    reserve = model[:reserve]
-
+    switch_off = model[:switch_off]
+    switch_on = model[:switch_on]
    T = model[:instance].time
    for t in 1:T
        # Startup limit
@@ -146,15 +95,8 @@ function _add_startup_shutdown_limit_eqs!(
        )
        # Shutdown limit
        if g.initial_power > g.shutdown_limit
-            # TODO check what happens with these variables when exporting the model
-            # Generator producing too much to be turned off in the first time period
-            # (can a binary variable have bounds x = 0?)
-            if formulation_status_vars.fix_vars_via_constraint
-                eq_shutdown_limit[g.name, 0] =
-                    @constraint(model, model[:switch_off][g.name, 1] <= 0.0)
-            else
-                fix(model[:switch_off][g.name, 1], 0.0; force = true)
-            end
+            eq_shutdown_limit[g.name, 0] =
+                @constraint(model, switch_off[g.name, 1] <= 0)
        end
        if t < T
            eq_shutdown_limit[g.name, t] = @constraint(
@@ -169,32 +111,6 @@ function _add_startup_shutdown_limit_eqs!(
    return
 end

-"""
-    _add_shutdown_cost_eqs!(model::JuMP.Model, g::Unit)::Nothing
-
-Variables
---
-* `switch_off`
-"""
-function _add_shutdown_cost_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    T = model[:instance].time
-    gi = g.name
-    for t in 1:T
-        shutdown_cost = 0.0
-        if shutdown_cost > 1e-7
-            # Equation (62) in Knueven et al. (2020)
-            add_to_expression!(
-                model[:obj],
-                model[:switch_off][gi, t],
-                shutdown_cost,
-            )
-        end
-    end # loop over time
-end
-
-"""
-    _add_ramp_eqs!(model, unit, formulation)
-"""
 function _add_ramp_eqs!(
    model::JuMP.Model,
    g::Unit,
@@ -241,24 +157,6 @@ function _add_ramp_eqs!(
    end
 end

-"""
-    _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
-
-Ensure constraints on up/down time are met.
-Based on Garver (1962), Malkin (2003), and Rajan and Takritti (2005).
-Eqns. (3), (4), (5) in Knueven et al. (2020).
-
-Variables
---
-* `is_on`
-* `switch_off`
-* `switch_on`
-
-Constraints
---
-* `eq_min_uptime`
-* `eq_min_downtime`
-"""
 function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
    is_on = model[:is_on]
    switch_off = model[:switch_off]
@@ -268,24 +166,18 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
    T = model[:instance].time
    for t in 1:T
        # Minimum up-time
-        # Equation (4) in Knueven et al. (2020)
        eq_min_uptime[g.name, t] = @constraint(
            model,
            sum(switch_on[g.name, i] for i in (t-g.min_uptime+1):t if i >= 1) <= is_on[g.name, t]
        )
-
        # Minimum down-time
-        # Equation (5) in Knueven et al. (2020)
        eq_min_downtime[g.name, t] = @constraint(
            model,
            sum(
                switch_off[g.name, i] for i in (t-g.min_downtime+1):t if i >= 1
            ) <= 1 - is_on[g.name, t]
        )
-
        # Minimum up/down-time for initial periods
-        # Equations (3a) and (3b) in Knueven et al. (2020)
-        # (using :switch_on and :switch_off instead of :is_on)
        if t == 1
            if g.initial_status > 0
                eq_min_uptime[g.name, 0] = @constraint(
@@ -308,9 +200,6 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
    end
 end

-"""
-    _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
-"""
 function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
    expr_net_injection = model[:expr_net_injection]
    for t in 1:model[:instance].time
--- a/src/transform/randomize.jl
+++ b/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!
--- a/src/transform/randomize/XavQiuAhm2021.jl
+++ b/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!
--- a/src/utils/benchmark.jl
+++ b/src/utils/benchmark.jl
@@ -0,0 +1,116 @@
+# 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 Distributed
+using Random
+
+function _run_benchmark_sample(;
+    case::String,
+    method::SolutionMethod,
+    formulation::Formulation,
+    solution_filename::String,
+    optimizer,
+)::Nothing
+    total_time = @elapsed begin
+        @info "Reading: $case"
+        time_read = @elapsed begin
+            instance = read_benchmark(case)
+        end
+        @info @sprintf("Read problem in %.2f seconds", time_read)
+        BLAS.set_num_threads(Threads.nthreads())
+        model = build_model(
+            instance = instance,
+            formulation = formulation,
+            optimizer = optimizer,
+            variable_names = true,
+        )
+        @info "Optimizing..."
+        BLAS.set_num_threads(1)
+        optimize!(model, method)
+    end
+    @info @sprintf("Total time was %.2f seconds", total_time)
+
+    @info "Writing solution: $solution_filename"
+    solution = UnitCommitment.solution(model)
+    write("$solution_filename", solution)
+
+    @info "Verifying solution..."
+    validate(instance, solution)
+    return
+end
+
+function _run_benchmark_combination(
+    case::String,
+    optimizer_name::String,
+    optimizer,
+    method_name::String,
+    method::SolutionMethod,
+    formulation_name::String,
+    formulation::Formulation,
+    trial,
+)
+    dirname = "results/$optimizer_name/$method_name/$formulation_name/$case"
+    function info(msg)
+        @info @sprintf(
+            "%-8s %-16s %-16s %-16s %-8s %s",
+            msg,
+            optimizer_name,
+            method_name,
+            formulation_name,
+            trial,
+            case
+        )
+    end
+    mkpath(dirname)
+    trial_filename = @sprintf("%s/%03d.json", dirname, trial)
+    if isfile(trial_filename)
+        info("skip")
+        return
+    end
+    info("run")
+    open("$trial_filename.log", "w") do file
+        redirect_stdout(file) do
+            redirect_stderr(file) do
+                return _run_benchmark_sample(
+                    case = case,
+                    method = method,
+                    formulation = formulation,
+                    solution_filename = trial_filename,
+                    optimizer = optimizer,
+                )
+            end
+        end
+    end
+    return info("done")
+end
+
+function _run_benchmarks(;
+    cases::Vector{String},
+    optimizers::Dict,
+    formulations::Dict,
+    methods::Dict,
+    trials,
+)
+    combinations = [
+        (c, s.first, s.second, m.first, m.second, f.first, f.second, t) for
+        c in cases for s in optimizers for f in formulations for
+        m in methods for t in trials
+    ]
+    shuffle!(combinations)
+    if nworkers() > 1
+        @printf("%24s", "")
+    end
+    @info @sprintf(
+        "%-8s %-16s %-16s %-16s %-8s %s",
+        "STATUS",
+        "SOLVER",
+        "METHOD",
+        "FORMULATION",
+        "TRIAL",
+        "CASE"
+    )
+    @sync @distributed for c in combinations
+        _run_benchmark_combination(c...)
+    end
+end
--- a/src/utils/sysimage.jl
+++ b/src/utils/sysimage.jl
@@ -1,28 +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 PackageCompiler
-
-using DataStructures
-using Distributions
-using JSON
-using JuMP
-using MathOptInterface
-using SparseArrays
-
-pkg = [
-    :DataStructures,
-    :Distributions,
-    :JSON,
-    :JuMP,
-    :MathOptInterface,
-    :SparseArrays,
-]
-
-@info "Building system image..."
-create_sysimage(
-    pkg,
-    precompile_statements_file = "build/precompile.jl",
-    sysimage_path = "build/sysimage.so",
-)
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -0,0 +1,26 @@
+[deps]
+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"
+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"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+DataStructures = "0.18"
+Distributions = "0.25"
+GZip = "0.5"
+JSON = "0.21"
+JuMP = "0.21"
+MathOptInterface = "0.9"
+PackageCompiler = "1"
+julia = "1"
--- a/test/import/egret_test.jl
+++ b/test/import/egret_test.jl
@@ -4,9 +4,12 @@

 using UnitCommitment

+basedir = @__DIR__
+
@testset "read_egret_solution" begin
-    solution =
-        UnitCommitment.read_egret_solution("fixtures/egret_output.json.gz")
+    solution = UnitCommitment.read_egret_solution(
+        "$basedir/../fixtures/egret_output.json.gz",
+    )
    for attr in ["Is on", "Production (MW)", "Production cost (\$)"]
        @test attr in keys(solution)
        @test "115_STEAM_1" in keys(solution[attr])
--- a/test/instance/read_test.jl
+++ b/test/instance/read_test.jl
@@ -22,6 +22,7 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
    @test instance.lines[5].normal_flow_limit == [1e8 for t in 1:4]
    @test instance.lines[5].emergency_flow_limit == [1e8 for t in 1:4]
    @test instance.lines[5].flow_limit_penalty == [5e3 for t in 1:4]
+    @test instance.lines_by_name["l5"].name == "l5"

    @test instance.lines[1].name == "l1"
    @test instance.lines[1].source.name == "b1"
@@ -34,6 +35,7 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip

    @test instance.buses[9].name == "b9"
    @test instance.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353]
+    @test instance.buses_by_name["b9"].name == "b9"

    unit = instance.units[1]
    @test unit.name == "g1"
@@ -62,6 +64,7 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
    @test unit.startup_categories[1].cost == 1000.0
    @test unit.startup_categories[2].cost == 1500.0
    @test unit.startup_categories[3].cost == 2000.0
+    @test instance.units_by_name["g1"].name == "g1"

    unit = instance.units[2]
    @test unit.name == "g2"
@@ -92,12 +95,15 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip

    @test instance.contingencies[1].lines == [instance.lines[1]]
    @test instance.contingencies[1].units == []
+    @test instance.contingencies[1].name == "c1"
+    @test instance.contingencies_by_name["c1"].name == "c1"

    load = instance.price_sensitive_loads[1]
    @test load.name == "ps1"
    @test load.bus.name == "b3"
    @test load.revenue == [100.0 for t in 1:4]
    @test load.demand == [50.0 for t in 1:4]
+    @test instance.price_sensitive_loads_by_name["ps1"].name == "ps1"
 end

@testset "read_benchmark sub-hourly" begin
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@
 using Test
 using UnitCommitment

+push!(Base.LOAD_PATH, @__DIR__)
 UnitCommitment._setup_logger()

 const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
@@ -28,7 +29,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")
--- a/test/transform/initcond_test.jl
+++ b/test/transform/initcond_test.jl
@@ -4,9 +4,12 @@

 using UnitCommitment, Cbc, JuMP

+basedir = @__DIR__
+
@testset "generate_initial_conditions!" begin
    # Load instance
-    instance = UnitCommitment.read("$(pwd())/fixtures/case118-initcond.json.gz")
+    instance =
+        UnitCommitment.read("$basedir/../fixtures/case118-initcond.json.gz")
    optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)

    # All units should have unknown initial conditions
--- a/test/transform/randomize/XavQiuAhm2021_test.jl
+++ b/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
--- a/test/transform/randomize_test.jl
+++ b/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
--- a/test/validation/repair_test.jl
+++ b/test/validation/repair_test.jl
@@ -4,9 +4,11 @@

 using UnitCommitment, JSON, GZip, DataStructures

+basedir = @__DIR__
+
 function parse_case14()
    return JSON.parse(
-        GZip.gzopen("../instances/test/case14.json.gz"),
+        GZip.gzopen("$basedir/../../instances/test/case14.json.gz"),
        dicttype = () -> DefaultOrderedDict(nothing),
    )
 end
Author	SHA1	Message	Date
Alinson S Xavier	93d3e5987d	Replace sysimage.jl by juliaw; add deps/formatter	2021-08-31 09:51:36 -05:00
Alinson S Xavier	f235333551	Improve benchmark scripts	2021-08-31 08:03:21 -05:00
Alinson S Xavier	6c566e0e79	Improve sysimage.jl	2021-08-20 04:51:09 -05:00
Alinson S Xavier	5c3f7b15d3	UnitCommitmentInstance: add _by_name fields	2021-08-19 07:07:25 -05:00
Alinson S Xavier	7c907a6eb5	Implement randomization method from XavQiuAhm2021	2021-08-05 17:04:37 -05:00
Alinson S. Xavier	b1498c50b3	GitHub Actions: Test fewer combinations	2021-07-26 07:57:17 -05:00
Aleksandr Kazachkov	000215e991	Add reserve shortfall penalty	2021-07-26 07:54:45 -05:00