From 3baddf158af31669eb545358ac5afe1ca3911254 Mon Sep 17 00:00:00 2001 From: Aleksandr Kazachkov Date: Fri, 20 Nov 2020 10:52:23 -0600 Subject: [PATCH] Implement multiple formulations --- Manifest.toml | 139 ++-- Project.toml | 1 + benchmark/Makefile | 241 +++++- benchmark/Manifest.toml | 170 ++-- benchmark/run.jl | 77 +- scripts/instances.txt | 182 +++++ scripts/run_batch.sh | 49 ++ src/UnitCommitment.jl | 7 +- src/components.jl | 46 ++ src/constraints.jl | 31 + src/formulation.jl | 1711 +++++++++++++++++++++++++++++++++++++++ src/instance.jl | 26 +- src/model2.jl | 475 +++++++++++ src/validate.jl | 27 +- src/variables.jl | 471 +++++++++++ test/model_test.jl | 62 +- 16 files changed, 3496 insertions(+), 219 deletions(-) create mode 100644 scripts/instances.txt create mode 100644 scripts/run_batch.sh create mode 100644 src/components.jl create mode 100644 src/constraints.jl create mode 100644 src/formulation.jl create mode 100644 src/model2.jl create mode 100644 src/variables.jl diff --git a/Manifest.toml b/Manifest.toml index c5b52a2..4368635 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,5 +1,11 @@ # This file is machine-generated - editing it directly is not advised +[[Artifacts]] +deps = ["Pkg"] +git-tree-sha1 = "c30985d8821e0cd73870b17b0ed0ce6dc44cb744" +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +version = "1.3.0" + [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" @@ -16,15 +22,15 @@ uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" version = "0.5.10" [[Bzip2_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "3663bfffede2ef41358b6fc2e1d8a6d50b3c3904" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e" uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" -version = "1.0.6+2" +version = "1.0.6+5" [[CEnum]] -git-tree-sha1 = "1b77a77c3b28e0b3f413f7567c9bb8dd9bdccd14" +git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9" uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" -version = "0.3.0" +version = "0.4.1" [[Calculus]] deps = ["LinearAlgebra"] @@ -34,9 +40,9 @@ version = "0.5.1" [[Cbc]] deps = ["BinaryProvider", "CEnum", "Cbc_jll", "Libdl", "MathOptInterface", "SparseArrays"] -git-tree-sha1 = "72e4299de0995a60a6230079adc7e47580870815" +git-tree-sha1 = "929d0500c50387e7ac7ae9956ca7d7ce5312c90d" uuid = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" -version = "0.7.0" +version = "0.7.1" [[Cbc_jll]] deps = ["Cgl_jll", "Clp_jll", "CoinUtils_jll", "CompilerSupportLibraries_jll", "Libdl", "OpenBLAS32_jll", "Osi_jll", "Pkg"] @@ -52,9 +58,9 @@ version = "0.60.2+5" [[Clp_jll]] deps = ["CoinUtils_jll", "CompilerSupportLibraries_jll", "Libdl", "OpenBLAS32_jll", "Osi_jll", "Pkg"] -git-tree-sha1 = "70fe9e52fd95fa37f645e3d30f08f436cc5b1457" +git-tree-sha1 = "79263d9383ca89b35f31c33ab5b880536a8413a4" uuid = "06985876-5285-5a41-9fcb-8948a742cc53" -version = "1.17.6+5" +version = "1.17.6+6" [[CodecBzip2]] deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"] @@ -80,22 +86,32 @@ git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" version = "0.3.0" +[[Compat]] +deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] +git-tree-sha1 = "48608f94f3e1a755f65e7ef34684675bb3653030" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "3.21.0" + [[CompilerSupportLibraries_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "8e695f735fca77e9708e795eda62afdb869cbb70" uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "0.3.3+0" +version = "0.3.4+0" [[DataStructures]] -deps = ["InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "edad9434967fdc0a2631a65d902228400642120c" +deps = ["Compat", "InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "db07bb22795762895b60e44d62b34b16c982a687" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.17.19" +version = "0.18.7" [[Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +[[DelimitedFiles]] +deps = ["Mmap"] +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" + [[DiffResults]] deps = ["StaticArrays"] git-tree-sha1 = "da24935df8e0c6cf28de340b958f6aac88eaa0cc" @@ -114,15 +130,15 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[DocStringExtensions]] deps = ["LibGit2", "Markdown", "Pkg", "Test"] -git-tree-sha1 = "c5714d9bcdba66389612dc4c47ed827c64112997" +git-tree-sha1 = "50ddf44c53698f5e784bbebb3f4b21c5807401b1" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.2" +version = "0.8.3" [[Documenter]] deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "1c593d1efa27437ed9dd365d1143c594b563e138" +git-tree-sha1 = "fb1ff838470573adc15c71ba79f8d31328f035da" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.25.1" +version = "0.25.2" [[ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"] @@ -131,8 +147,8 @@ uuid = "f6369f11-7733-5829-9624-2563aa707210" version = "0.10.12" [[GLPK]] -deps = ["BinaryProvider", "GLPK_jll", "Libdl", "MathOptInterface", "SparseArrays"] -git-tree-sha1 = "86573ecb852e303b209212046af44871f1c0e49c" +deps = ["BinaryProvider", "Libdl", "MathOptInterface", "SparseArrays"] +git-tree-sha1 = "3420033e843e140d9237238d69937a5bc7292e5a" uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" version = "0.13.0" @@ -143,10 +159,10 @@ uuid = "e8aa6df9-e6ca-548a-97ff-1f85fc5b8b98" version = "4.64.0+0" [[GMP_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "4dd9301d3a027c05ec403e756ee7a60e3c367e5d" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "15abc5f976569a1c9d651aff02f7222ef305eb2a" uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d" -version = "6.1.2+5" +version = "6.1.2+6" [[GZip]] deps = ["Libdl"] @@ -156,9 +172,9 @@ version = "0.5.1" [[HTTP]] deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"] -git-tree-sha1 = "2ac03263ce44be4222342bca1c51c36ce7566161" +git-tree-sha1 = "c7ec02c4c6a039a98a15f955462cd7aea5df4508" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.8.17" +version = "0.8.19" [[IniFile]] deps = ["Test"] @@ -170,23 +186,28 @@ version = "0.5.0" deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[JLLWrappers]] +git-tree-sha1 = "7cec881362e5b4e367ff0279dd99a06526d51a55" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.1.2" + [[JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e" +git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.0" +version = "0.21.1" [[JSONSchema]] deps = ["HTTP", "JSON", "ZipFile"] -git-tree-sha1 = "832a4d327d9dafdae55a6ecae04f9997c83615cc" +git-tree-sha1 = "a9ecdbc90be216912a2e3e8a8e38dc4c93f0d065" uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" -version = "0.3.0" +version = "0.3.2" [[JuMP]] -deps = ["Calculus", "DataStructures", "ForwardDiff", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Random", "SparseArrays", "Statistics"] -git-tree-sha1 = "cbab42e2e912109d27046aa88f02a283a9abac7c" +deps = ["Calculus", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Random", "SparseArrays", "Statistics"] +git-tree-sha1 = "57c17a221a55f81890aabf00f478886859e25eaf" uuid = "4076af6c-e467-56ae-b986-b466b2749572" -version = "0.21.3" +version = "0.21.5" [[LibGit2]] deps = ["Printf"] @@ -204,9 +225,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[MacroTools]] deps = ["Markdown", "Random"] -git-tree-sha1 = "f7d2e3f654af75f01ec49be82c231c382214223a" +git-tree-sha1 = "6a8a2a625ab0dea913aba95c11370589e0239ff0" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.5" +version = "0.5.6" [[Markdown]] deps = ["Base64"] @@ -214,9 +235,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" [[MathOptInterface]] deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"] -git-tree-sha1 = "cd2049c055c7d192a235670d50faa375361624ba" +git-tree-sha1 = "5a1d631e0a9087d425e024d66b9c71e92e78fda8" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -version = "0.9.14" +version = "0.9.17" [[MbedTLS]] deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"] @@ -226,9 +247,9 @@ version = "1.0.2" [[MbedTLS_jll]] deps = ["Libdl", "Pkg"] -git-tree-sha1 = "a0cb0d489819fa7ea5f9fa84c7e7eba19d8073af" +git-tree-sha1 = "c0b1286883cac4e2b617539de41111e0776d02e8" uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.16.6+1" +version = "2.16.8+0" [[Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" @@ -245,21 +266,21 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "0.3.4" [[OpenBLAS32_jll]] -deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"] -git-tree-sha1 = "793b33911239d2651c356c823492b58d6490d36a" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19c33675cdeb572c1b17f96c492459d4f4958036" uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2" -version = "0.3.9+4" +version = "0.3.10+0" [[OpenSpecFun_jll]] -deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"] -git-tree-sha1 = "d51c416559217d974a1113522d5919235ae67a87" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "9db77584158d0ab52307f8c04f8e7c08ca76b5b3" uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.3+3" +version = "0.5.3+4" [[OrderedCollections]] -git-tree-sha1 = "293b70ac1780f9584c89268a6e2a560d938a7065" +git-tree-sha1 = "16c08bf5dba06609fe45e30860092d6fa41fde7b" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.3.0" +version = "1.3.1" [[Osi_jll]] deps = ["CoinUtils_jll", "CompilerSupportLibraries_jll", "Libdl", "OpenBLAS32_jll", "Pkg"] @@ -274,10 +295,10 @@ uuid = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" version = "1.2.1" [[Parsers]] -deps = ["Dates", "Test"] -git-tree-sha1 = "10134f2ee0b1978ae7752c41306e131a684e1f06" +deps = ["Dates"] +git-tree-sha1 = "6fa4202675c05ba0f8268a6ddf07606350eda3ce" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.0.7" +version = "1.0.11" [[Pkg]] deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] @@ -297,9 +318,9 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[Requires]] deps = ["UUIDs"] -git-tree-sha1 = "d37400976e98018ee840e0ca4f9d20baa231dc6b" +git-tree-sha1 = "28faf1c963ca1dc3ec87f166d92982e3c4a1f66d" uuid = "ae029012-a4dd-5104-9daa-d747884805df" -version = "1.0.1" +version = "1.1.0" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" @@ -307,6 +328,10 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" [[Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +[[SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + [[Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" @@ -355,12 +380,12 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" [[ZipFile]] deps = ["Libdl", "Printf", "Zlib_jll"] -git-tree-sha1 = "254975fef2fc526583bb9b7c9420fe66ffe09f2f" +git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7" uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" -version = "0.9.2" +version = "0.9.3" [[Zlib_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "320228915c8debb12cb434c59057290f0834dbf6" uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.11+14" +version = "1.2.11+18" diff --git a/Project.toml b/Project.toml index 272768d..5daf318 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Requires = "ae029012-a4dd-5104-9daa-d747884805df" diff --git a/benchmark/Makefile b/benchmark/Makefile index be4c39b..4665bea 100644 --- a/benchmark/Makefile +++ b/benchmark/Makefile @@ -6,6 +6,9 @@ SHELL := /bin/bash JULIA := julia --project=. --sysimage ../build/sysimage.so TIMESTAMP := $(shell date "+%Y-%m-%d %H:%M") SRC_FILES := $(wildcard ../src/*.jl) +DEST := . +FORMULATION := tight +results_dir := results_$(FORMULATION) INSTANCES_PGLIB := \ pglib-uc/ca/2014-09-01_reserves_0 \ @@ -87,11 +90,205 @@ INSTANCES_MATPOWER := \ matpower/case6515rte/2017-02-01 \ matpower/case6515rte/2017-08-01 -SAMPLES := 1 2 3 -SOLUTIONS_MATPOWER := $(foreach s,$(SAMPLES),$(addprefix results/,$(addsuffix .$(s).sol.json,$(INSTANCES_MATPOWER)))) -SOLUTIONS_PGLIB := $(foreach s,$(SAMPLES),$(addprefix results/,$(addsuffix .$(s).sol.json,$(INSTANCES_PGLIB)))) +INSTANCES_INFORMS1 := \ + matpower/case1888rte/2017-01-01 \ + matpower/case1888rte/2017-01-02 \ + matpower/case1888rte/2017-01-03 \ + matpower/case1888rte/2017-01-04 \ + matpower/case1888rte/2017-01-05 \ + matpower/case1888rte/2017-01-06 \ + matpower/case1888rte/2017-01-07 \ + matpower/case1888rte/2017-01-08 \ + matpower/case1888rte/2017-01-09 \ + matpower/case1888rte/2017-01-10 \ + matpower/case1888rte/2017-01-11 \ + matpower/case1888rte/2017-01-12 \ + matpower/case1888rte/2017-01-13 \ + matpower/case1888rte/2017-01-14 \ + matpower/case1888rte/2017-01-15 \ + matpower/case1888rte/2017-01-16 \ + matpower/case1888rte/2017-01-17 \ + matpower/case1888rte/2017-01-18 \ + matpower/case1888rte/2017-01-19 \ + matpower/case1888rte/2017-01-20 \ + matpower/case1888rte/2017-01-21 \ + matpower/case1888rte/2017-01-22 \ + matpower/case1888rte/2017-01-23 \ + matpower/case1888rte/2017-01-24 \ + matpower/case1888rte/2017-01-25 \ + matpower/case1888rte/2017-01-26 \ + matpower/case1888rte/2017-01-27 \ + matpower/case1888rte/2017-01-28 \ + matpower/case1888rte/2017-01-29 \ + matpower/case1888rte/2017-01-30 \ + matpower/case1888rte/2017-01-31 \ + matpower/case1888rte/2017-02-01 \ + matpower/case1888rte/2017-02-02 \ + matpower/case1888rte/2017-02-03 \ + matpower/case1888rte/2017-02-04 \ + matpower/case1888rte/2017-02-05 \ + matpower/case1888rte/2017-02-06 \ + matpower/case1888rte/2017-02-07 \ + matpower/case1888rte/2017-02-08 \ + matpower/case1888rte/2017-02-09 \ + matpower/case1888rte/2017-02-10 \ + matpower/case1888rte/2017-02-11 \ + matpower/case1888rte/2017-02-12 \ + matpower/case1888rte/2017-02-13 \ + matpower/case1888rte/2017-02-14 \ + matpower/case1888rte/2017-02-15 \ + matpower/case1888rte/2017-02-16 \ + matpower/case1888rte/2017-02-17 \ + matpower/case1888rte/2017-02-18 \ + matpower/case1888rte/2017-02-19 \ + matpower/case1888rte/2017-02-20 \ + matpower/case1888rte/2017-02-21 \ + matpower/case1888rte/2017-02-22 \ + matpower/case1888rte/2017-02-23 \ + matpower/case1888rte/2017-02-24 \ + matpower/case1888rte/2017-02-25 \ + matpower/case1888rte/2017-02-26 \ + matpower/case1888rte/2017-02-27 \ + matpower/case1888rte/2017-02-28 \ + matpower/case1888rte/2017-03-01 -.PHONY: tables save small large clean-mps matpower pglib +INSTANCES_INFORMS2 := \ + matpower/case3375wp/2017-01-01 \ + matpower/case3375wp/2017-01-02 \ + matpower/case3375wp/2017-01-03 \ + matpower/case3375wp/2017-01-04 \ + matpower/case3375wp/2017-01-05 \ + matpower/case3375wp/2017-01-06 \ + matpower/case3375wp/2017-01-07 \ + matpower/case3375wp/2017-01-08 \ + matpower/case3375wp/2017-01-09 \ + matpower/case3375wp/2017-01-10 \ + matpower/case3375wp/2017-01-11 \ + matpower/case3375wp/2017-01-12 \ + matpower/case3375wp/2017-01-13 \ + matpower/case3375wp/2017-01-14 \ + matpower/case3375wp/2017-01-15 \ + matpower/case3375wp/2017-01-16 \ + matpower/case3375wp/2017-01-17 \ + matpower/case3375wp/2017-01-18 \ + matpower/case3375wp/2017-01-19 \ + matpower/case3375wp/2017-01-20 \ + matpower/case3375wp/2017-01-21 \ + matpower/case3375wp/2017-01-22 \ + matpower/case3375wp/2017-01-23 \ + matpower/case3375wp/2017-01-24 \ + matpower/case3375wp/2017-01-25 \ + matpower/case3375wp/2017-01-26 \ + matpower/case3375wp/2017-01-27 \ + matpower/case3375wp/2017-01-28 \ + matpower/case3375wp/2017-01-29 \ + matpower/case3375wp/2017-01-30 \ + matpower/case3375wp/2017-01-31 \ + matpower/case3375wp/2017-02-01 \ + matpower/case3375wp/2017-02-02 \ + matpower/case3375wp/2017-02-03 \ + matpower/case3375wp/2017-02-04 \ + matpower/case3375wp/2017-02-05 \ + matpower/case3375wp/2017-02-06 \ + matpower/case3375wp/2017-02-07 \ + matpower/case3375wp/2017-02-08 \ + matpower/case3375wp/2017-02-09 \ + matpower/case3375wp/2017-02-10 \ + matpower/case3375wp/2017-02-11 \ + matpower/case3375wp/2017-02-12 \ + matpower/case3375wp/2017-02-13 \ + matpower/case3375wp/2017-02-14 \ + matpower/case3375wp/2017-02-15 \ + matpower/case3375wp/2017-02-16 \ + matpower/case3375wp/2017-02-17 \ + matpower/case3375wp/2017-02-18 \ + matpower/case3375wp/2017-02-19 \ + matpower/case3375wp/2017-02-20 \ + matpower/case3375wp/2017-02-21 \ + matpower/case3375wp/2017-02-22 \ + matpower/case3375wp/2017-02-23 \ + matpower/case3375wp/2017-02-24 \ + matpower/case3375wp/2017-02-25 \ + matpower/case3375wp/2017-02-26 \ + matpower/case3375wp/2017-02-27 \ + matpower/case3375wp/2017-02-28 \ + matpower/case3375wp/2017-03-01 + +INSTANCES_INFORMS3 := \ + matpower/case6468rte/2017-01-01 \ + matpower/case6468rte/2017-01-02 \ + matpower/case6468rte/2017-01-03 \ + matpower/case6468rte/2017-01-04 \ + matpower/case6468rte/2017-01-05 \ + matpower/case6468rte/2017-01-06 \ + matpower/case6468rte/2017-01-07 \ + matpower/case6468rte/2017-01-08 \ + matpower/case6468rte/2017-01-09 \ + matpower/case6468rte/2017-01-10 \ + matpower/case6468rte/2017-01-11 \ + matpower/case6468rte/2017-01-12 \ + matpower/case6468rte/2017-01-13 \ + matpower/case6468rte/2017-01-14 \ + matpower/case6468rte/2017-01-15 \ + matpower/case6468rte/2017-01-16 \ + matpower/case6468rte/2017-01-17 \ + matpower/case6468rte/2017-01-18 \ + matpower/case6468rte/2017-01-19 \ + matpower/case6468rte/2017-01-20 \ + matpower/case6468rte/2017-01-21 \ + matpower/case6468rte/2017-01-22 \ + matpower/case6468rte/2017-01-23 \ + matpower/case6468rte/2017-01-24 \ + matpower/case6468rte/2017-01-25 \ + matpower/case6468rte/2017-01-26 \ + matpower/case6468rte/2017-01-27 \ + matpower/case6468rte/2017-01-28 \ + matpower/case6468rte/2017-01-29 \ + matpower/case6468rte/2017-01-30 \ + matpower/case6468rte/2017-01-31 \ + matpower/case6468rte/2017-02-01 \ + matpower/case6468rte/2017-02-02 \ + matpower/case6468rte/2017-02-03 \ + matpower/case6468rte/2017-02-04 \ + matpower/case6468rte/2017-02-05 \ + matpower/case6468rte/2017-02-06 \ + matpower/case6468rte/2017-02-07 \ + matpower/case6468rte/2017-02-08 \ + matpower/case6468rte/2017-02-09 \ + matpower/case6468rte/2017-02-10 \ + matpower/case6468rte/2017-02-11 \ + matpower/case6468rte/2017-02-12 \ + matpower/case6468rte/2017-02-13 \ + matpower/case6468rte/2017-02-14 \ + matpower/case6468rte/2017-02-15 \ + matpower/case6468rte/2017-02-16 \ + matpower/case6468rte/2017-02-17 \ + matpower/case6468rte/2017-02-18 \ + matpower/case6468rte/2017-02-19 \ + matpower/case6468rte/2017-02-20 \ + matpower/case6468rte/2017-02-21 \ + matpower/case6468rte/2017-02-22 \ + matpower/case6468rte/2017-02-23 \ + matpower/case6468rte/2017-02-24 \ + matpower/case6468rte/2017-02-25 \ + matpower/case6468rte/2017-02-26 \ + matpower/case6468rte/2017-02-27 \ + matpower/case6468rte/2017-02-28 \ + matpower/case6468rte/2017-03-01 + +INSTANCES_TEST := \ + test/case14 + +#SAMPLES := 1 2 3 +SAMPLES := 1 +SOLUTIONS_MATPOWER := $(foreach s,$(SAMPLES),$(addprefix $(results_dir)/,$(addsuffix .$(s).sol.json,$(INSTANCES_MATPOWER)))) +SOLUTIONS_PGLIB := $(foreach s,$(SAMPLES),$(addprefix $(results_dir)/,$(addsuffix .$(s).sol.json,$(INSTANCES_PGLIB)))) +SOLUTIONS_INFORMS1 := $(foreach s,$(SAMPLES),$(addprefix $(results_dir)/,$(addsuffix .$(s).sol.json,$(INSTANCES_INFORMS1)))) +SOLUTIONS_INFORMS2 := $(foreach s,$(SAMPLES),$(addprefix $(results_dir)/,$(addsuffix .$(s).sol.json,$(INSTANCES_INFORMS2)))) +SOLUTIONS_INFORMS3 := $(foreach s,$(SAMPLES),$(addprefix $(results_dir)/,$(addsuffix .$(s).sol.json,$(INSTANCES_INFORMS3)))) +SOLUTIONS_TEST := $(foreach s,$(SAMPLES),$(addprefix $(results_dir)/,$(addsuffix .$(s).sol.json,$(INSTANCES_TEST)))) + +.PHONY: tables save small large clean-mps matpower pglib informs1 informs2 informs3 test all: matpower pglib @@ -99,23 +296,47 @@ matpower: $(SOLUTIONS_MATPOWER) pglib: $(SOLUTIONS_PGLIB) +informs1: $(SOLUTIONS_INFORMS1) +informs2: $(SOLUTIONS_INFORMS2) +informs3: $(SOLUTIONS_INFORMS3) + +test: $(SOLUTIONS_TEST) + clean: - @rm -rf tables/benchmark* tables/compare* results + @rm -rf tables/benchmark* tables/compare* $(results_dir) clean-mps: - @rm -fv results/*/*/*.mps.gz + @rm -fv $(results_dir)/*/*/*.mps.gz clean-sol: - @rm -rf results/*/*/*.sol.* + @rm -rf $(results_dir)/*/*/*.sol.* save: mkdir -p "runs/$(TIMESTAMP)" - rsync -avP results tables "runs/$(TIMESTAMP)/" + rsync -avP $(results_dir) tables "runs/$(TIMESTAMP)/" results/%.sol.json: run.jl @echo "run $*" - @mkdir -p $(dir results/$*) - @$(JULIA) run.jl $* 2>&1 | cat > results/$*.log + @mkdir -p $(dir $(DEST)/$(results_dir)/$*) + @$(JULIA) run.jl $* default $(DEST)/$(results_dir) 2>&1 | cat > $(DEST)/$(results_dir)/$*.log + @echo "run $* [done]" + +results_tight/%.sol.json: run.jl + @echo "run $*" + @mkdir -p $(dir $(DEST)/$(results_dir)/$*) + @$(JULIA) run.jl $* tight $(DEST)/$(results_dir) 2>&1 | cat > $(DEST)/$(results_dir)/$*.log + @echo "run $* [done]" + +results_default/%.sol.json: run.jl + @echo "run $*" + @mkdir -p $(dir $(DEST)/$(results_dir)/$*) + @$(JULIA) run.jl $* default $(DEST)/$(results_dir) 2>&1 | cat > $(DEST)/$(results_dir)/$*.log + @echo "run $* [done]" + +results_sparse/%.sol.json: run.jl + @echo "run $*" + @mkdir -p $(dir $(DEST)/$(results_dir)/$*) + @$(JULIA) run.jl $* sparse $(DEST)/$(results_dir) 2>&1 | cat > $(DEST)/$(results_dir)/$*.log @echo "run $* [done]" tables: diff --git a/benchmark/Manifest.toml b/benchmark/Manifest.toml index c319e25..40357a7 100644 --- a/benchmark/Manifest.toml +++ b/benchmark/Manifest.toml @@ -1,5 +1,11 @@ # This file is machine-generated - editing it directly is not advised +[[Artifacts]] +deps = ["Pkg"] +git-tree-sha1 = "c30985d8821e0cd73870b17b0ed0ce6dc44cb744" +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +version = "1.3.0" + [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" @@ -16,15 +22,15 @@ uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" version = "0.5.10" [[Bzip2_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "3663bfffede2ef41358b6fc2e1d8a6d50b3c3904" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e" uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" -version = "1.0.6+2" +version = "1.0.6+5" [[CEnum]] -git-tree-sha1 = "1b77a77c3b28e0b3f413f7567c9bb8dd9bdccd14" +git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9" uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" -version = "0.3.0" +version = "0.4.1" [[Calculus]] deps = ["LinearAlgebra"] @@ -34,9 +40,9 @@ version = "0.5.1" [[Cbc]] deps = ["BinaryProvider", "CEnum", "Cbc_jll", "Libdl", "MathOptInterface", "SparseArrays"] -git-tree-sha1 = "72e4299de0995a60a6230079adc7e47580870815" +git-tree-sha1 = "929d0500c50387e7ac7ae9956ca7d7ce5312c90d" uuid = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" -version = "0.7.0" +version = "0.7.1" [[Cbc_jll]] deps = ["Cgl_jll", "Clp_jll", "CoinUtils_jll", "CompilerSupportLibraries_jll", "Libdl", "OpenBLAS32_jll", "Osi_jll", "Pkg"] @@ -52,15 +58,9 @@ version = "0.60.2+5" [[Clp_jll]] deps = ["CoinUtils_jll", "CompilerSupportLibraries_jll", "Libdl", "OpenBLAS32_jll", "Osi_jll", "Pkg"] -git-tree-sha1 = "70fe9e52fd95fa37f645e3d30f08f436cc5b1457" +git-tree-sha1 = "79263d9383ca89b35f31c33ab5b880536a8413a4" uuid = "06985876-5285-5a41-9fcb-8948a742cc53" -version = "1.17.6+5" - -[[CodeTracking]] -deps = ["InteractiveUtils", "UUIDs"] -git-tree-sha1 = "cab4da992adc0a64f63fa30d2db2fd8bec40cab4" -uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" -version = "0.5.11" +version = "1.17.6+6" [[CodecBzip2]] deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"] @@ -87,16 +87,16 @@ uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" version = "0.3.0" [[CompilerSupportLibraries_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "8e695f735fca77e9708e795eda62afdb869cbb70" uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "0.3.3+0" +version = "0.3.4+0" [[DataStructures]] deps = ["InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "edad9434967fdc0a2631a65d902228400642120c" +git-tree-sha1 = "88d48e133e6d3dd68183309877eac74393daa7eb" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.17.19" +version = "0.17.20" [[Dates]] deps = ["Printf"] @@ -120,18 +120,15 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[DocStringExtensions]] deps = ["LibGit2", "Markdown", "Pkg", "Test"] -git-tree-sha1 = "c5714d9bcdba66389612dc4c47ed827c64112997" +git-tree-sha1 = "50ddf44c53698f5e784bbebb3f4b21c5807401b1" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.2" +version = "0.8.3" [[Documenter]] deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "1c593d1efa27437ed9dd365d1143c594b563e138" +git-tree-sha1 = "fb1ff838470573adc15c71ba79f8d31328f035da" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.25.1" - -[[FileWatching]] -uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" +version = "0.25.2" [[ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"] @@ -140,10 +137,10 @@ uuid = "f6369f11-7733-5829-9624-2563aa707210" version = "0.10.12" [[GLPK]] -deps = ["BinaryProvider", "GLPK_jll", "Libdl", "MathOptInterface", "SparseArrays"] -git-tree-sha1 = "86573ecb852e303b209212046af44871f1c0e49c" +deps = ["BinaryProvider", "CEnum", "GLPK_jll", "Libdl", "MathOptInterface"] +git-tree-sha1 = "0984f1669480cdecd465458b4abf81b238fbfe50" uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" -version = "0.13.0" +version = "0.14.2" [[GLPK_jll]] deps = ["GMP_jll", "Libdl", "Pkg"] @@ -152,10 +149,10 @@ uuid = "e8aa6df9-e6ca-548a-97ff-1f85fc5b8b98" version = "4.64.0+0" [[GMP_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "4dd9301d3a027c05ec403e756ee7a60e3c367e5d" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "15abc5f976569a1c9d651aff02f7222ef305eb2a" uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d" -version = "6.1.2+5" +version = "6.1.2+6" [[GZip]] deps = ["Libdl"] @@ -164,16 +161,16 @@ uuid = "92fee26a-97fe-5a0c-ad85-20a5f3185b63" version = "0.5.1" [[Gurobi]] -deps = ["Libdl", "LinearAlgebra", "MathOptInterface", "MathProgBase", "SparseArrays"] -git-tree-sha1 = "f36a2fa62909675681aec582ccfc4a4a629406e4" +deps = ["CEnum", "Libdl", "MathOptInterface"] +git-tree-sha1 = "de2015da3bffcf005ef51b78163e81bfb7b2301d" uuid = "2e9cd046-0924-5485-92f1-d5272153d98b" -version = "0.8.1" +version = "0.9.2" [[HTTP]] deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"] -git-tree-sha1 = "2ac03263ce44be4222342bca1c51c36ce7566161" +git-tree-sha1 = "c7ec02c4c6a039a98a15f955462cd7aea5df4508" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.8.17" +version = "0.8.19" [[IniFile]] deps = ["Test"] @@ -185,6 +182,11 @@ version = "0.5.0" deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[JLLWrappers]] +git-tree-sha1 = "c70593677bbf2c3ccab4f7500d0f4dacfff7b75c" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.1.3" + [[JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e" @@ -193,9 +195,9 @@ version = "0.21.0" [[JSONSchema]] deps = ["HTTP", "JSON", "ZipFile"] -git-tree-sha1 = "832a4d327d9dafdae55a6ecae04f9997c83615cc" +git-tree-sha1 = "a9ecdbc90be216912a2e3e8a8e38dc4c93f0d065" uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" -version = "0.3.0" +version = "0.3.2" [[JuMP]] deps = ["Calculus", "DataStructures", "ForwardDiff", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Random", "SparseArrays", "Statistics"] @@ -203,12 +205,6 @@ git-tree-sha1 = "cbab42e2e912109d27046aa88f02a283a9abac7c" uuid = "4076af6c-e467-56ae-b986-b466b2749572" version = "0.21.3" -[[JuliaInterpreter]] -deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"] -git-tree-sha1 = "79e4496b79e8af45198f8c291f26d4514d6b06d6" -uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a" -version = "0.7.24" - [[LibGit2]] deps = ["Printf"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" @@ -223,51 +219,33 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" -[[LoweredCodeUtils]] -deps = ["JuliaInterpreter"] -git-tree-sha1 = "1b632dc108106101a9909db7be8f8b32ed8d02f7" -uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b" -version = "0.4.6" - [[MacroTools]] deps = ["Markdown", "Random"] -git-tree-sha1 = "f7d2e3f654af75f01ec49be82c231c382214223a" +git-tree-sha1 = "6a8a2a625ab0dea913aba95c11370589e0239ff0" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.5" +version = "0.5.6" [[Markdown]] deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" -[[MathOptFormat]] -deps = ["CodecZlib", "DataStructures", "HTTP", "JSON", "JSONSchema", "MathOptInterface"] -git-tree-sha1 = "0206edd9310b863c222af23f71fde5d09e95c20c" -uuid = "f4570300-c277-12e8-125c-4912f86ce65d" -version = "0.2.2" - [[MathOptInterface]] deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"] -git-tree-sha1 = "cd2049c055c7d192a235670d50faa375361624ba" +git-tree-sha1 = "5a1d631e0a9087d425e024d66b9c71e92e78fda8" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -version = "0.9.14" - -[[MathProgBase]] -deps = ["LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "9abbe463a1e9fc507f12a69e7f29346c2cdc472c" -uuid = "fdba3010-5040-5b88-9595-932c9decdf73" -version = "0.7.8" +version = "0.9.17" [[MbedTLS]] deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"] -git-tree-sha1 = "426a6978b03a97ceb7ead77775a1da066343ec6e" +git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe" uuid = "739be429-bea8-5141-9913-cc70e7f3736d" -version = "1.0.2" +version = "1.0.3" [[MbedTLS_jll]] deps = ["Libdl", "Pkg"] -git-tree-sha1 = "a0cb0d489819fa7ea5f9fa84c7e7eba19d8073af" +git-tree-sha1 = "c0b1286883cac4e2b617539de41111e0776d02e8" uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.16.6+1" +version = "2.16.8+0" [[Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" @@ -284,21 +262,21 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "0.3.4" [[OpenBLAS32_jll]] -deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"] -git-tree-sha1 = "793b33911239d2651c356c823492b58d6490d36a" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19c33675cdeb572c1b17f96c492459d4f4958036" uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2" -version = "0.3.9+4" +version = "0.3.10+0" [[OpenSpecFun_jll]] -deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"] -git-tree-sha1 = "d51c416559217d974a1113522d5919235ae67a87" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "9db77584158d0ab52307f8c04f8e7c08ca76b5b3" uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.3+3" +version = "0.5.3+4" [[OrderedCollections]] -git-tree-sha1 = "293b70ac1780f9584c89268a6e2a560d938a7065" +git-tree-sha1 = "16c08bf5dba06609fe45e30860092d6fa41fde7b" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.3.0" +version = "1.3.1" [[Osi_jll]] deps = ["CoinUtils_jll", "CompilerSupportLibraries_jll", "Libdl", "OpenBLAS32_jll", "Pkg"] @@ -308,15 +286,15 @@ version = "0.108.5+3" [[PackageCompiler]] deps = ["Libdl", "Pkg", "UUIDs"] -git-tree-sha1 = "98aa9c653e1dc3473bb5050caf8501293db9eee1" +git-tree-sha1 = "3eee77c94646163f15bd8626acf494360897f890" uuid = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" -version = "1.2.1" +version = "1.2.3" [[Parsers]] -deps = ["Dates", "Test"] -git-tree-sha1 = "10134f2ee0b1978ae7752c41306e131a684e1f06" +deps = ["Dates"] +git-tree-sha1 = "6fa4202675c05ba0f8268a6ddf07606350eda3ce" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.0.7" +version = "1.0.11" [[Pkg]] deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] @@ -336,15 +314,9 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[Requires]] deps = ["UUIDs"] -git-tree-sha1 = "d37400976e98018ee840e0ca4f9d20baa231dc6b" +git-tree-sha1 = "28faf1c963ca1dc3ec87f166d92982e3c4a1f66d" uuid = "ae029012-a4dd-5104-9daa-d747884805df" -version = "1.0.1" - -[[Revise]] -deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "UUIDs", "Unicode"] -git-tree-sha1 = "0992d4643e27b2deb9f2e4ec7a56b7033813a027" -uuid = "295af30f-e4ad-537b-8983-00126c2a3abe" -version = "2.7.3" +version = "1.1.0" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" @@ -399,19 +371,19 @@ uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" [[UnitCommitment]] -deps = ["Cbc", "DataStructures", "Documenter", "GLPK", "GZip", "JSON", "JuMP", "LinearAlgebra", "Logging", "MathOptFormat", "MathOptInterface", "PackageCompiler", "Printf", "Requires", "Revise", "SparseArrays", "Test", "TimerOutputs"] +deps = ["Cbc", "DataStructures", "Documenter", "GLPK", "GZip", "Gurobi", "JSON", "JuMP", "LinearAlgebra", "Logging", "MathOptInterface", "OrderedCollections", "PackageCompiler", "Printf", "Requires", "SparseArrays", "Test", "TimerOutputs"] path = ".." uuid = "64606440-39ea-11e9-0f29-3303a1d3d877" version = "2.1.0" [[ZipFile]] deps = ["Libdl", "Printf", "Zlib_jll"] -git-tree-sha1 = "254975fef2fc526583bb9b7c9420fe66ffe09f2f" +git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7" uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" -version = "0.9.2" +version = "0.9.3" [[Zlib_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "320228915c8debb12cb434c59057290f0834dbf6" uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.11+14" +version = "1.2.11+18" diff --git a/benchmark/run.jl b/benchmark/run.jl index b3c6fb9..a757e30 100644 --- a/benchmark/run.jl +++ b/benchmark/run.jl @@ -11,14 +11,47 @@ using Printf using LinearAlgebra function main() - basename, suffix = split(ARGS[1], ".") - solution_filename = "results/$basename.$suffix.sol.json" - model_filename = "results/$basename.$suffix.mps.gz" - + NUM_THREADS = 4 time_limit = 60 * 20 + BLAS.set_num_threads(NUM_THREADS) + + if length(ARGS) >= 2 + mode = string("_", ARGS[2]) + else + mode = "_default" + end + if length(ARGS) >= 3 && !isempty(strip(ARGS[3])) + results_dir = ARGS[3] + else + results_dir = string("./","results$mode") + end + + # Validate mode and set formulation + if mode == "_default" + formulation = UnitCommitment.DefaultFormulation + elseif mode == "_tight" + formulation = UnitCommitment.TightFormulation + elseif mode == "_sparse" + formulation = UnitCommitment.SparseDefaultFormulation + else + error("Unknown formulation requested: ", ARGS[2]) + end + + # Filename is instance_name.sample_number.sol.gz + # Parse out the instance + sample parts to create output files + basename, suffix = split(ARGS[1], ".") # will not work if suffix part is not present + model_filename_stub = string(results_dir,"/$basename.$suffix") + solution_filename = string("$model_filename_stub.sol.json") - BLAS.set_num_threads(4) - global_logger(TimeLogger(initial_time = time())) + # Choose logging options + logname, logfile = nothing, nothing + #logname = string("$model_filename_stub.out") + if isa(logname, String) && !isempty(logname) + logfile = open(logname, "w") + global_logger(TimeLogger(initial_time = time(), file = logfile)) + else + global_logger(TimeLogger(initial_time = time())) + end total_time = @elapsed begin @info "Reading: $basename" @@ -28,18 +61,26 @@ function main() @info @sprintf("Read problem in %.2f seconds", time_read) time_model = @elapsed begin - model = build_model(instance=instance, - optimizer=optimizer_with_attributes(Gurobi.Optimizer, - "Threads" => 4, - "Seed" => rand(1:1000), - )) + optimizer=optimizer_with_attributes(Gurobi.Optimizer, + "Threads" => NUM_THREADS, + "Seed" => rand(1:1000)) + model = build_model(instance=instance, optimizer=optimizer, formulation=formulation) end + end + @info "Setting names..." + UnitCommitment.set_variable_names!(model) + + model_filename = string(model_filename_stub,".init",".mps.gz") + @info string("Exporting initial model without transmission constraints to ", model_filename) + JuMP.write_to_file(model.mip, model_filename) + + total_time += @elapsed begin @info "Optimizing..." BLAS.set_num_threads(1) UnitCommitment.optimize!(model, time_limit=time_limit, gap_limit=1e-3) - end + @info @sprintf("Total time was %.2f seconds", total_time) @info "Writing: $solution_filename" @@ -51,11 +92,13 @@ function main() @info "Verifying solution..." UnitCommitment.validate(instance, solution) - @info "Setting variable names..." - UnitCommitment.set_variable_names!(model) - - @info "Exporting model..." + model_filename = string(model_filename_stub,".final",".mps.gz") + @info string("Exporting final model to ", model_filename) JuMP.write_to_file(model.mip, model_filename) -end + + if !isnothing(logfile) + close(logfile) + end +end # main main() diff --git a/scripts/instances.txt b/scripts/instances.txt new file mode 100644 index 0000000..3527baf --- /dev/null +++ b/scripts/instances.txt @@ -0,0 +1,182 @@ +matpower/case1888rte/2017-01-01 +matpower/case1888rte/2017-01-02 +matpower/case1888rte/2017-01-03 +matpower/case1888rte/2017-01-04 +matpower/case1888rte/2017-01-05 +matpower/case1888rte/2017-01-06 +matpower/case1888rte/2017-01-07 +matpower/case1888rte/2017-01-08 +matpower/case1888rte/2017-01-09 +matpower/case1888rte/2017-01-10 +matpower/case1888rte/2017-01-11 +matpower/case1888rte/2017-01-12 +matpower/case1888rte/2017-01-13 +matpower/case1888rte/2017-01-14 +matpower/case1888rte/2017-01-15 +matpower/case1888rte/2017-01-16 +matpower/case1888rte/2017-01-17 +matpower/case1888rte/2017-01-18 +matpower/case1888rte/2017-01-19 +matpower/case1888rte/2017-01-20 +matpower/case1888rte/2017-01-21 +matpower/case1888rte/2017-01-22 +matpower/case1888rte/2017-01-23 +matpower/case1888rte/2017-01-24 +matpower/case1888rte/2017-01-25 +matpower/case1888rte/2017-01-26 +matpower/case1888rte/2017-01-27 +matpower/case1888rte/2017-01-28 +matpower/case1888rte/2017-01-29 +matpower/case1888rte/2017-01-30 +matpower/case1888rte/2017-01-31 +matpower/case1888rte/2017-02-01 +matpower/case1888rte/2017-02-02 +matpower/case1888rte/2017-02-03 +matpower/case1888rte/2017-02-04 +matpower/case1888rte/2017-02-05 +matpower/case1888rte/2017-02-06 +matpower/case1888rte/2017-02-07 +matpower/case1888rte/2017-02-08 +matpower/case1888rte/2017-02-09 +matpower/case1888rte/2017-02-10 +matpower/case1888rte/2017-02-11 +matpower/case1888rte/2017-02-12 +matpower/case1888rte/2017-02-13 +matpower/case1888rte/2017-02-14 +matpower/case1888rte/2017-02-15 +matpower/case1888rte/2017-02-16 +matpower/case1888rte/2017-02-17 +matpower/case1888rte/2017-02-18 +matpower/case1888rte/2017-02-19 +matpower/case1888rte/2017-02-20 +matpower/case1888rte/2017-02-21 +matpower/case1888rte/2017-02-22 +matpower/case1888rte/2017-02-23 +matpower/case1888rte/2017-02-24 +matpower/case1888rte/2017-02-25 +matpower/case1888rte/2017-02-26 +matpower/case1888rte/2017-02-27 +matpower/case1888rte/2017-02-28 +matpower/case1888rte/2017-03-01 +matpower/case3375wp/2017-01-01 +matpower/case3375wp/2017-01-02 +matpower/case3375wp/2017-01-03 +matpower/case3375wp/2017-01-04 +matpower/case3375wp/2017-01-05 +matpower/case3375wp/2017-01-06 +matpower/case3375wp/2017-01-07 +matpower/case3375wp/2017-01-08 +matpower/case3375wp/2017-01-09 +matpower/case3375wp/2017-01-10 +matpower/case3375wp/2017-01-11 +matpower/case3375wp/2017-01-12 +matpower/case3375wp/2017-01-13 +matpower/case3375wp/2017-01-14 +matpower/case3375wp/2017-01-15 +matpower/case3375wp/2017-01-16 +matpower/case3375wp/2017-01-17 +matpower/case3375wp/2017-01-18 +matpower/case3375wp/2017-01-19 +matpower/case3375wp/2017-01-20 +matpower/case3375wp/2017-01-21 +matpower/case3375wp/2017-01-22 +matpower/case3375wp/2017-01-23 +matpower/case3375wp/2017-01-24 +matpower/case3375wp/2017-01-25 +matpower/case3375wp/2017-01-26 +matpower/case3375wp/2017-01-27 +matpower/case3375wp/2017-01-28 +matpower/case3375wp/2017-01-29 +matpower/case3375wp/2017-01-30 +matpower/case3375wp/2017-01-31 +matpower/case3375wp/2017-02-01 +matpower/case3375wp/2017-02-02 +matpower/case3375wp/2017-02-03 +matpower/case3375wp/2017-02-04 +matpower/case3375wp/2017-02-05 +matpower/case3375wp/2017-02-06 +matpower/case3375wp/2017-02-07 +matpower/case3375wp/2017-02-08 +matpower/case3375wp/2017-02-09 +matpower/case3375wp/2017-02-10 +matpower/case3375wp/2017-02-11 +matpower/case3375wp/2017-02-12 +matpower/case3375wp/2017-02-13 +matpower/case3375wp/2017-02-14 +matpower/case3375wp/2017-02-15 +matpower/case3375wp/2017-02-16 +matpower/case3375wp/2017-02-17 +matpower/case3375wp/2017-02-18 +matpower/case3375wp/2017-02-19 +matpower/case3375wp/2017-02-20 +matpower/case3375wp/2017-02-21 +matpower/case3375wp/2017-02-22 +matpower/case3375wp/2017-02-23 +matpower/case3375wp/2017-02-24 +matpower/case3375wp/2017-02-25 +matpower/case3375wp/2017-02-26 +matpower/case3375wp/2017-02-27 +matpower/case3375wp/2017-02-28 +matpower/case3375wp/2017-03-01 +matpower/case6468rte/2017-01-01 +matpower/case6468rte/2017-01-02 +matpower/case6468rte/2017-01-03 +matpower/case6468rte/2017-01-04 +matpower/case6468rte/2017-01-05 +matpower/case6468rte/2017-01-06 +matpower/case6468rte/2017-01-07 +matpower/case6468rte/2017-01-08 +matpower/case6468rte/2017-01-09 +matpower/case6468rte/2017-01-10 +matpower/case6468rte/2017-01-11 +matpower/case6468rte/2017-01-12 +matpower/case6468rte/2017-01-13 +matpower/case6468rte/2017-01-14 +matpower/case6468rte/2017-01-15 +matpower/case6468rte/2017-01-16 +matpower/case6468rte/2017-01-17 +matpower/case6468rte/2017-01-18 +matpower/case6468rte/2017-01-19 +matpower/case6468rte/2017-01-20 +matpower/case6468rte/2017-01-21 +matpower/case6468rte/2017-01-22 +matpower/case6468rte/2017-01-23 +matpower/case6468rte/2017-01-24 +matpower/case6468rte/2017-01-25 +matpower/case6468rte/2017-01-26 +matpower/case6468rte/2017-01-27 +matpower/case6468rte/2017-01-28 +matpower/case6468rte/2017-01-29 +matpower/case6468rte/2017-01-30 +matpower/case6468rte/2017-01-31 +matpower/case6468rte/2017-02-01 +matpower/case6468rte/2017-02-02 +matpower/case6468rte/2017-02-03 +matpower/case6468rte/2017-02-04 +matpower/case6468rte/2017-02-05 +matpower/case6468rte/2017-02-06 +matpower/case6468rte/2017-02-07 +matpower/case6468rte/2017-02-08 +matpower/case6468rte/2017-02-09 +matpower/case6468rte/2017-02-10 +matpower/case6468rte/2017-02-11 +matpower/case6468rte/2017-02-12 +matpower/case6468rte/2017-02-13 +matpower/case6468rte/2017-02-14 +matpower/case6468rte/2017-02-15 +matpower/case6468rte/2017-02-16 +matpower/case6468rte/2017-02-17 +matpower/case6468rte/2017-02-18 +matpower/case6468rte/2017-02-19 +matpower/case6468rte/2017-02-20 +matpower/case6468rte/2017-02-21 +matpower/case6468rte/2017-02-22 +matpower/case6468rte/2017-02-23 +matpower/case6468rte/2017-02-24 +matpower/case6468rte/2017-02-25 +matpower/case6468rte/2017-02-26 +matpower/case6468rte/2017-02-27 +matpower/case6468rte/2017-02-28 +matpower/case6468rte/2017-03-01 + +test/case14 diff --git a/scripts/run_batch.sh b/scripts/run_batch.sh new file mode 100644 index 0000000..54a4cd8 --- /dev/null +++ b/scripts/run_batch.sh @@ -0,0 +1,49 @@ +#!/bin/bash +#SBATCH --array=1-180 +#SBATCH --time=02:00:00 +#SBATCH --account=def-alodi +#SBATCH --mem-per-cpu=1G +#SBATCH --cpus-per-task=4 +#SBATCH --mail-user=aleksandr.kazachkov@polymtl.ca +#SBATCH --mail-type=BEGIN +#SBATCH --mail-type=END +#SBATCH --mail-type=FAIL +#SBATCH --array=182 +#SBATCH --time=00:00:30 +#SBATCH --mem-per-cpu=500M +#SBATCH --cpus-per-task=1 +#SBATCH --time=01:00:00 +#SBATCH --mem-per-cpu=1G +#SBATCH --cpus-per-task=4 + +MODE="tight" +if [ ! -z $1 ]; then + MODE=$1 +fi + +#CASE_NUM=`printf %03d $SLURM_ARRAY_TASK_ID` +PROJ_DIR="${REPOS_DIR}/UnitCommitment2.jl" +INST=$(sed -n "${SLURM_ARRAY_TASK_ID}p" ${PROJ_DIR}/scripts/instances.txt) +#DEST="${PROJ_DIR}/benchmark" +DEST="${HOME}/scratch/uc" +RESULTS_DIR="${DEST}/results_${MODE}" +NUM_SAMPLES=1 + +if [ $MODE == "sparse" ] || [ $MODE == "default" ] || [ $MODE == "tight" ] +then + echo "Running task $SLURM_ARRAY_TASK_ID for instance $INST with results sent to ${RESULTS_DIR}" +else + echo "Unrecognized mode: $1. Exiting." + exit +fi + +cd ${PROJ_DIR}/benchmark +mkdir -p $(dirname ${RESULTS_DIR}/${INST}) +for i in $(seq ${NUM_SAMPLES}); do + FILE=$INST.$i + #echo "Running $FILE at `date` using command julia --project=${PROJ_DIR}/benchmark --sysimage=${PROJ_DIR}/build/sysimage.so ${PROJ_DIR}/benchmark/run.jl ${FILE} ${MODE} ${RESULTS_DIR} 2&>1 | cat > ${RESULTS_DIR}/${FILE}.log" + #julia --project=${PROJ_DIR}/benchmark --sysimage=${PROJ_DIR}/build/sysimage.so ${PROJ_DIR}/benchmark/run.jl ${FILE} ${MODE} ${RESULTS_DIR} 2&>1 | cat > ${RESULTS_DIR}/${FILE}.log + echo "Running $FILE at `date` using command julia --project=${PROJ_DIR}/benchmark --sysimage=${PROJ_DIR}/build/sysimage.so ${PROJ_DIR}/benchmark/run.jl ${FILE} ${MODE} ${RESULTS_DIR} &> ${RESULTS_DIR}/${FILE}.log" + julia --project=${PROJ_DIR}/benchmark --sysimage=${PROJ_DIR}/build/sysimage.so ${PROJ_DIR}/benchmark/run.jl ${FILE} ${MODE} ${RESULTS_DIR} &> ${RESULTS_DIR}/${FILE}.log + #julia --project=${PROJ_DIR}/benchmark --sysimage=${PROJ_DIR}/build/sysimage.so ${PROJ_DIR}/benchmark/run.jl ${FILE} ${MODE} ${RESULTS_DIR} +done diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 7c335d3..d37e790 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -7,7 +7,12 @@ module UnitCommitment include("dotdict.jl") include("instance.jl") include("screening.jl") - include("model.jl") + include("components.jl") + include("variables.jl") + include("constraints.jl") + include("formulation.jl") + #include("model.jl") + include("model2.jl") include("sensitivity.jl") include("validate.jl") include("convert.jl") diff --git a/src/components.jl b/src/components.jl new file mode 100644 index 0000000..7f7f2a4 --- /dev/null +++ b/src/components.jl @@ -0,0 +1,46 @@ +################################################## +# Component types +abstract type UCComponentType end +abstract type RequiredConstraints <: UCComponentType end +abstract type SystemConstraints <: UCComponentType end +abstract type GenerationLimits <: UCComponentType end +abstract type PiecewiseProduction <: UCComponentType end +abstract type UpDownTime <: UCComponentType end +abstract type ReserveConstraints <: UCComponentType end +abstract type RampLimits <: UCComponentType end +abstract type StartupCosts <: UCComponentType end +abstract type ShutdownCosts <: UCComponentType end + +################################################## +# Components +""" +Generic component of the unit commitment problem. + +Elements +=== + * `name`: name of the component + * `description`: gives a brief summary of what the component adds + * `type`: reference back to the UCComponentType being modeled + * `vars`: required variables + * `constrs`: constraints that are created by this function + * `add_component`: function to add constraints and update the objective to capture this component + * `params`: extra parameters the component might use +""" +mutable struct UCComponent + "Name of the component." + name::String + "Description of what the component adds." + description::String + "Which part of the unit commitment problem is modeled by this component." + type::Type{<:UCComponentType} + "Variables that are needed for the component (subset of `var_list`)." + vars::Union{Array{Symbol},Nothing} + "Equations that are modified for the component (subset of `constr_list`)." + constrs::Union{Array{Symbol},Nothing} + "Function to add constraints and objective coefficients needed for this component to the model. Signature should be (component, mip, model)." + add_component::Union{Function,Nothing} + "Extra parameters for the component." + params::Any +end # struct UCComponent + +export UCComponent diff --git a/src/constraints.jl b/src/constraints.jl new file mode 100644 index 0000000..5612052 --- /dev/null +++ b/src/constraints.jl @@ -0,0 +1,31 @@ +################################################## +# Constraints +""" +List of constraints that the model will potentially have +""" +constr_list = + [ + :startup_choose, + :startup_restrict, + :segprod_limit, + :segprod_limita, + :segprod_limitb, + :prod_above_def, + :prod_limit, + :str_prod_limit, + :binary_link, + :switch_on_off, + :ramp_up, + :ramp_down, + :str_ramp_up, + :str_ramp_down, + :startstop_limit, + :startup_limit, + :shutdown_limit, + :min_uptime, + :min_downtime, + :power_balance, + :net_injection_def, + :min_reserve + ] + diff --git a/src/formulation.jl b/src/formulation.jl new file mode 100644 index 0000000..bbbf2a0 --- /dev/null +++ b/src/formulation.jl @@ -0,0 +1,1711 @@ +# Determine which formulation is being used + +using DataStructures # for OrderedDict +using JuMP + +RESERVES_WHEN_START_UP = true +RESERVES_WHEN_RAMP_UP = true +RESERVES_WHEN_RAMP_DOWN = true +RESERVES_WHEN_SHUT_DOWN = true + +mutable struct UnitCommitmentModel2 + mip::JuMP.Model + vars::DotDict + eqs::DotDict + exprs::DotDict + instance::UnitCommitmentInstance + isf::Array{Float64, 2} + lodf::Array{Float64, 2} + obj::AffExpr + components::Array{UCComponent,1} +end # UnitCommitmentModel2 + + +################################################## +## Required constraints + +function add_required_constraints_default(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + T = model.instance.time + for g in model.instance.units + gi = g.name + if !all(g.must_run) && any(g.must_run) + error("Partially must-run units are not currently supported.") + end + + known_initial_conditions = (g.initial_status != nothing && g.initial_power != nothing) + if known_initial_conditions + if g.initial_status > 0 && g.initial_power > 0 + # 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(model.vars.switch_on[gi, 1], 0.0; force = true) + + # amk added below + # (redundant to min_uptime[gi, 0] constraint) + # Unit will continue to be on for min_uptime - g.initial_status periods + UT = max(g.min_uptime - g.initial_status, 0) + for t = 1:min(UT, T) + fix(model.vars.switch_off[gi, t], 0.0; force = true) + fix(model.vars.is_on[gi, t], 1.0; force = true) + if t < T + # If on in t, then cannot be turned on in t+1 + fix(model.vars.switch_on[gi, t+1], 0.0; force = true) + end + end + end + + if g.initial_status < 0 + # Generator is initially off (for -g.initial_status time periods) + # Cannot be switched off more + fix(model.vars.switch_off[gi, 1], 0.0; force = true) + + # amk added below, but probably redundant due to min_downtime[gi, 0] constraint + DT = max(g.min_downtime + g.initial_status, 0) + for t = 1:min(DT, T) + fix(model.vars.switch_on[gi, t], 0.0; force = true) + fix(model.vars.is_on[gi, t], 0.0; force = true) + if t < T + # If off in t, then cannot be turned off in t+1 + fix(model.vars.switch_off[gi, t+1], 0.0; force = true) + end + end + end + + if g.initial_power > g.shutdown_limit + # Generator producing too much to be turned off in the first time period + # TODO check what happens with these variables when exporting the model + # (can a binary variable have bounds x = 0?) + #eqs.shutdown_limit[gi, 0] = @constraint(mip, vars.switch_off[gi, 1] <= 0) + fix(vars.switch_off[gi, 1], 0.; force = true) + end + end # known_initial_conditions + + for t = 1:T + if !g.provides_spinning_reserves[t] + fix(model.vars.reserve[gi, t], 0.0; force = true) + end + + if g.must_run[t] + # If the generator _must_ run, then it is obviously on and cannot be switched off + fix(model.vars.is_on[gi, t], 1.0; force = true) + fix(model.vars.switch_off[gi, t], 0.0; force = true) + if t == 1 && known_initial_conditions && g.initial_status < 0 + # In the first time period, force to switch on if was off before + fix(model.vars.switch_on[gi, t], 1.0; force = true) + else + # Otherwise, it is on, and will never turn off, so will never need to turn on + fix(model.vars.switch_on[gi, t], 0.0; force = true) + end + end # g.must_run[t] + end # loop over time + end # loop over units + + for t = 1:T + for lm in model.instance.lines + add_to_expression!(model.obj, model.vars.overflow[lm.name, t], lm.flow_limit_penalty[t]) + end # loop over lines + end # loop over time +end # add_required_constraints_default + +DefaultRequiredConstraints = UCComponent( + "DefaultRequiredConstraints", + "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.", + RequiredConstraints, + [:is_on, :switch_on, :switch_off, :reserve, :overflow], + nothing, + add_required_constraints_default, + nothing +) # DefaultRequiredConstraints + + +################################################## +## System constraints + +function add_system_constraints_default(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + #nop +end # add_system_constraints_default + +DefaultSystemConstraints = UCComponent( + "DefaultSystemConstraints", + "Ensure constraints on system, across units / buses, are met.", + SystemConstraints, + nothing, + nothing, + add_system_constraints_default, + nothing +) # DefaultSystemConstraints + + +################################################## +## Generation limits + +function add_generation_limits_default(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for t = 1:T + for g in model.instance.units + gi = g.name + + # Objective function terms for production costs + # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term + add_to_expression!(model.obj, vars.is_on[gi, t], g.min_power_cost[t]) + + # Production limit + # Equation (18) in Kneuven et al. (2020) + # as \bar{p}_g(t) \le \bar{P}_g u_g(t) + # amk: this is a weaker version of (20) and (21) in Kneuven et al. (2020) + # but keeping it here in case those are not present + power_diff = max(g.max_power[t], 0.) - max(g.min_power[t], 0.) + if power_diff < 1e-7 + power_diff = 0. + end + eqs.prod_limit[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + vars.reserve[gi, t] + <= power_diff * vars.is_on[gi, t]) + end # loop over units + end # loop over time +end # add_generation_limits_default + +# Default +DefaultGenerationLimits = UCComponent( + "DefaultSystemConstraints", + "Ensure constraints on system are met."* + " Based on Garver (1962) and Morales-España et al. (2013)."* + " Eqns. (18), part of (69) in Kneuven et al. (2020).", + GenerationLimits, + [:is_on, :prod_above, :reserve], + [:prod_limit], + add_generation_limits_default, + nothing +) # DefaultGenerationLimits + + +################################################## +## Piecewise production + +function add_piecewise_production_Garver62(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for t = 1:T + for g in model.instance.units + gi = g.name + K = length(g.cost_segments) + for k in 1:K + # 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 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* + eqs.segprod_limit[gi, k, t] = + @constraint(mip, + vars.segprod[gi, k, t] + <= g.cost_segments[k].mw[t] * vars.is_on[gi, t]) + + # Also add this as an explicit upper bound on segprod to make the solver's work a bit easier + set_upper_bound(vars.segprod[gi, k, t], g.cost_segments[k].mw[t]) + + # Definition of production + # Equation (43) in Kneuven et al. (2020) + eqs.prod_above_def[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + == sum(vars.segprod[gi, k, t] for k in 1:K)) + + # Objective function + # Equation (44) in Kneuven et al. (2020) + add_to_expression!(model.obj, vars.segprod[gi, k, t], g.cost_segments[k].cost[t]) + end # loop over cost segments + end # loop over units + end # loop over time +end # add_piecewise_production_Garver62 + +# Garver, 1962 +PiecewiseProduction_Garver62 = UCComponent( + "PiecewiseProduction_Garver62", + "Ensure respect of production limits along each segment."* + " Based on Garver (1962)."* + " Equations (42), (43), (44) 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 that segment*.", + PiecewiseProduction, + [:segprod, :is_on, :prod_above], + [:segprod_limit, :prod_above_def], + add_piecewise_production_Garver62, + nothing +) # PiecewiseProduction_Garver62 + + +function add_piecewise_production_CarArr06(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for t = 1:T + for g in model.instance.units + gi = g.name + K = length(g.cost_segments) + for k in 1:K + # 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 that segment* + eqs.segprod_limit[gi, k, t] = + @constraint(mip, + vars.segprod[gi, k, t] + <= g.cost_segments[k].mw[t]) + + # Also add this as an explicit upper bound on segprod to make the solver's work a bit easier + set_upper_bound(vars.segprod[gi, k, t], g.cost_segments[k].mw[t]) + + # Definition of production + # Equation (43) in Kneuven et al. (2020) + eqs.prod_above_def[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + == sum(vars.segprod[gi, k, t] for k in 1:K)) + + # Objective function + # Equation (44) in Kneuven et al. (2020) + add_to_expression!(model.obj, vars.segprod[gi, k, t], g.cost_segments[k].cost[t]) + end # loop over cost segments + end # loop over units + end # loop over time +end # add_piecewise_production_CarArr06 + +# Carrion and Arroyo (2006) +PiecewiseProduction_CarArr06 = UCComponent( + "PiecewiseProduction_CarArr06", + "Ensure respect of production limits along each segment."* + " Based on Garver (1962) and Carrión and Arryo (2006),"* + " which replaces (42) in Kneuven et al. (2020) with a weaker version missing the on/off variable."* + " Equations (45), (43), (44) 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 that segment*.", + PiecewiseProduction, + [:segprod, :is_on, :prod_above], + [:segprod_limit, :prod_above_def], + add_piecewise_production_CarArr06, + nothing +) # PiecewiseProduction_CarArr06 + + +function add_piecewise_production_KneOstWat18(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for t = 1:T + for g in model.instance.units + gi = g.name + K = length(g.cost_segments) + for k in 1:K + # Pbar^{k-1) + Pbar0 = g.min_power[t] + ( k > 1 ? sum( g.cost_segments[ell].mw[t] for ell in 1:k-1 ) : 0. ) + # Pbar^k + Pbar1 = g.cost_segments[k].mw[t] + Pbar0 + + Cv = 0. + SU = g.startup_limit # startup rate + if Pbar1 <= SU + Cv = 0. + elseif Pbar0 < SU # && Pbar1 > SU + Cv = Pbar1 - SU + else # Pbar0 >= SU + Cv = g.cost_segments[k].mw[t] # this will imply that we cannot produce along this segment if switch_on = 1 + end + + Cw = 0. + SD = g.shutdown_limit # shutdown rate + if Pbar1 <= SD + Cw = 0. + elseif Pbar0 < SD # && Pbar1 > SD + Cw = Pbar1 - SD + else # Pbar0 >= SD + Cw = g.cost_segments[k].mw[t] + end + + + if g.min_uptime > 1 + # Equation (46) in Kneuven et al. (2020) + eqs.segprod_limit[gi, k, t] = + @constraint(mip, + vars.segprod[gi, k, t] + <= g.cost_segments[k].mw[t] * vars.is_on[gi, t] + - Cv * vars.switch_on[gi, t] + - (t < T ? Cw * vars.switch_off[gi, t+1] : 0.) + ) + else + # Equation (47a)/(48a) in Kneuven et al. (2020) + eqs.segprod_limita[gi, k, t] = + @constraint(mip, + vars.segprod[gi, k, t] + <= g.cost_segments[k].mw[t] * vars.is_on[gi, t] + - Cv * vars.switch_on[gi, t] + - (t < T ? max(0, Cv-Cw) * vars.switch_off[gi, t+1] : 0.) + ) + + # Equation (47b)/(48b) in Kneuven et al. (2020) + eqs.segprod_limitb[gi, k, t] = + @constraint(mip, + vars.segprod[gi, k, t] + <= g.cost_segments[k].mw[t] * vars.is_on[gi, t] + - max(0, Cw-Cv) * vars.switch_on[gi, t] + - (t < T ? Cw * vars.switch_off[gi, t+1] : 0.) + ) + end # check if g.min_uptime > 1 + + # Definition of production + # Equation (43) in Kneuven et al. (2020) + eqs.prod_above_def[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + == sum(vars.segprod[gi, k, t] for k in 1:K)) + + # Objective function + # Equation (44) in Kneuven et al. (2020) + add_to_expression!(model.obj, vars.segprod[gi, k, t], g.cost_segments[k].cost[t]) + + # Also add an explicit upper bound on segprod to make the solver's work a bit easier + set_upper_bound(vars.segprod[gi, k, t], g.cost_segments[k].mw[t]) + end # loop over cost segments + end # loop over units + end # loop over time +end # add_piecewise_production_KneOstWat18 + +# Replace (42) with (46) and (48) +PiecewiseProduction_KneOstWat18 = UCComponent( + "PiecewiseProduction_KneOstWat18", + "Ensure respect of production limits along each segment."* + " Based on Kneuven et al. (2018b)."* + " Eqns. (43), (44), (46), (48) 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 that segment*.", + PiecewiseProduction, + [:segprod, :is_on, :prod_above], + [:segprod_limit, :segprod_limita, :segprod_limitb, :prod_above_def], + add_piecewise_production_KneOstWat18, + nothing +) # PiecewiseProduction_KneOstWat18 + + + +################################################## +## Up / down time + +function add_updowntime_default(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for t = 1:T + for g in model.instance.units + if g.must_run[t] + continue + end + + gi = g.name + known_initial_conditions = (g.initial_status != nothing && g.initial_power != nothing) + is_initially_on = (known_initial_conditions && (g.initial_status > 0)) ? 1.0 : 0.0 + + # Link binary variables + # Equation (2) in Kneuven et al. (2020), originally from Garver (1962) + if t == 1 + if known_initial_conditions + # When initial conditions are unknown, we allow the generator to have is_on = 1 without switch_on = 1 + # This means we do not pay the startup cost in the first time period + eqs.binary_link[gi, t] = + @constraint(mip, + vars.is_on[gi, t] - is_initially_on == + vars.switch_on[gi, t] - vars.switch_off[gi, t]) + end + else + eqs.binary_link[gi, t] = + @constraint(mip, + vars.is_on[gi, t] - vars.is_on[gi, t-1] == + vars.switch_on[gi, t] - vars.switch_off[gi, t]) + end + + # Cannot switch on and off at the same time + # amk: I am not sure this is in Kneuven et al. (2020) + eqs.switch_on_off[gi, t] = + @constraint(mip, + vars.switch_on[gi, t] + vars.switch_off[gi, t] <= 1) + + # Minimum up/down-time for initial periods + # Equations (3a) and (3b) in Kneuven et al. (2020) + # (using :switch_on and :switch_off instead of :is_on) + if t == 1 && known_initial_conditions + if g.initial_status > 0 + eqs.min_uptime[gi, 0] = + @constraint(mip, sum(vars.switch_off[gi, i] + for i in 1:(g.min_uptime - g.initial_status) if i <= T) == 0) + else + eqs.min_downtime[gi, 0] = + @constraint(mip, sum(vars.switch_on[gi, i] + for i in 1:(g.min_downtime + g.initial_status) if i <= T) == 0) + end + end + + # Minimum up-time + # Equation (4) in Kneuven et al. (2020) + eqs.min_uptime[gi, t] = + @constraint(mip, + sum(vars.switch_on[gi, i] + for i in (t - g.min_uptime + 1):t if i >= 1 + ) <= vars.is_on[gi, t]) + + # Minimum down-time + # Equation (5) in Kneuven et al. (2020) + eqs.min_downtime[gi, t] = + @constraint(mip, + sum(vars.switch_off[gi, i] + for i in (t - g.min_downtime + 1):t if i >= 1 + ) <= 1 - vars.is_on[gi, t]) + end # loop over units + end # loop over time +end # add_updowntime_default + +# Default +DefaultUpDownTime = UCComponent( + "DefaultUpDownTime", + "Ensure constraints on up/down time are met."* + " Based on Garver (1962), Malkin (2003), and Rajan and Takritti (2005)."* + " Eqns. (2), (4), (5) in Kneuven et al. (2020).", + UpDownTime, + [:is_on, :switch_on, :switch_off, :reserve_shortfall], + [:binary_link, :min_uptime, :min_downtime, :switch_on_off], + add_updowntime_default, + nothing +) # DefaultUpDownTime + + +################################################## +## Reserves + +function add_reserves_default(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + T = model.instance.time + for t = 1:T + for g in model.instance.units + gi = g.name + if !g.provides_spinning_reserves[t] + fix(model.vars.reserve[gi, t], 0.0; force=true) + #else + #add_to_expression!(mip.exprs.reserve[g.bus.name, t], + # model.vars.reserve[gi, t], 1.0) + end + end # loop over units + + # Equation (68) in Kneuven et al. (2020) + # As in Morales-España et al. (2013a) + # Akin to the alternative formulation with max_power_avail + # from Carrión and Arroyo (2006) and Ostrowski et al. (2012) + shortfall_penalty = model.instance.shortfall_penalty[t] + model.eqs.min_reserve[t] = + @constraint(model.mip, + sum(model.vars.reserve[g.name, t] for g in model.instance.units) + + (shortfall_penalty > 1e-7 ? model.vars.reserve_shortfall[t] : 0.) # amk added + >= model.instance.reserves.spinning[t]) + + # amk added: Account for shortfall contribution to objective + if shortfall_penalty > 1e-7 + add_to_expression!(mip.obj, + shortfall_penalty, + model.vars.reserve_shortfall[t]) + else + # Not added to the model at all + #fix(model.vars.reserve_shortfall[t], 0.; force=true) + end + end # loop over time +end # add_reserves_default + +# Default reserves formulation +DefaultReserves = UCComponent( + "DefaultReserves", + "Ensure constraints on reserves are met."* + " Based on Morales-España et al. (2013a)."* + " Eqn. (68) from Kneuven et al. (2020).", + ReserveConstraints, + [:reserve, :reserve_shortfall], + [:min_reserve], + add_reserves_default, + nothing +) # DefaultReserves + + +function add_reserves_max_power_avail(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + error("add_reserves_max_power_avail is not implemented.") + + T = model.instance.time + for t = 1:T + for g in model.instance.units + gi = g.name + if !g.provides_spinning_reserves[t] + end + end # loop over units + + # Equation (67) in Kneuven et al. (2020) + # From Carrión and Arroyo (2006) and Ostrowski et al. (2012) + model.eqs.min_reserve[t] = + @constraint(model.mip, + model.vars.max_power_avail[gi, t] + >= sum(bus.load for bus in model.instance.buses) + model.instance.reserves.spinning[t] + ) + end # loop over time +end # add_reserves_max_power_avail + +# TODO not finished +# Reserves formulation using :max_power_avail +ReservesMaxPowerAvail = UCComponent( + "ReservesMaxPowerAvail", + "Ensure constraints on reserves are met."* + " Based on Carrión and Arroyo (2006) and Ostrowski et al. (2012)."* + " Eqn. (67) from Kneuven et al. (2020).", + ReserveConstraints, + [:max_power_avail, :reserve_shortfall], + [:min_reserve], + add_reserves_max_power_avail, + nothing +) # ReservesMaxPowerAvail + + +################################################## +## Ramping limits + +function add_ramping_ArrCon00(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + T = model.instance.time + vars, eqs = model.vars, model.eqs + for t = 1:T + for g in model.instance.units + gi = g.name + known_initial_conditions = (g.initial_status != nothing && g.initial_power != nothing) + is_initially_on = known_initial_conditions && (g.initial_status > 0) + RU = g.ramp_up_limit + RD = g.ramp_down_limit + SU = g.startup_limit + SD = g.shutdown_limit + + # Ramp up limit + if t == 1 + # Ignore ramping limits in first period if initial conditions are unknown + if known_initial_conditions && is_initially_on + # min power is _not_ multiplied by is_on because if !is_on, then ramp up is irrelevant + eqs.ramp_up[gi, t] = + @constraint(mip, + g.min_power[t] + vars.prod_above[gi, t] + + (RESERVES_WHEN_RAMP_UP ? vars.reserve[gi, t] : 0.) + <= g.initial_power + RU) + end + else + max_prod_this_period = g.min_power[t] * vars.is_on[gi, t] + vars.prod_above[gi, t] + max_prod_this_period += ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? vars.reserve[gi, t] : 0. ) + min_prod_last_period = g.min_power[t-1] * vars.is_on[gi, t-1] + vars.prod_above[gi, t-1] + + # Equation (24) in Kneuven et al. (2020) + eqs.ramp_up[gi, t] = + @constraint(mip, + max_prod_this_period - min_prod_last_period + <= RU * vars.is_on[gi, t-1] + SU * vars.switch_on[gi, t]) + end # check if t = 1 or not + + # Ramp down limit + if t == 1 + # Ignore ramping limits in first period if initial conditions are unknown + if known_initial_conditions && is_initially_on + # TODO If RD < SD, or more specifically if + # min_power + RD < initial_power < SD + # then the generator should be able to shut down at time t = 1, + # but the constraint below will force the unit to produce power + eqs.ramp_down[gi, t] = + @constraint(mip, + g.initial_power + - (g.min_power[t] + vars.prod_above[gi, t]) + <= RD) + end + else + max_prod_last_period = g.min_power[t-1] * vars.is_on[gi, t-1] + vars.prod_above[gi, t-1] + max_prod_last_period += ( RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? vars.reserve[gi, t-1] : 0. ) # amk added + min_prod_this_period = g.min_power[t] * vars.is_on[gi, t] + vars.prod_above[gi, t] + + # Equation (25) in Kneuven et al. (2020) + eqs.ramp_down[gi, t] = + @constraint(mip, + max_prod_last_period - min_prod_this_period + <= RD * vars.is_on[gi, t] + SD * vars.switch_off[gi, t]) + end # check if t = 1 or not + end # loop over units + end # loop over time +end # add_ramping_ArrCon00 + +# Denser version of MorEsp13 ramping +Ramping_ArrCon00 = UCComponent( + "Ramping_ArrCon00", + "Ensure constraints on ramping are met."* + " Based on Arroyo and Conejo (2000)."* + " Eqns. (24), (25) in Kneuven et al. (2020).", + RampLimits, + [:is_on, :prod_above, :reserve], + [:ramp_up, :ramp_down], + add_ramping_ArrCon00, + nothing +) # Ramping_ArrCon00 + + +function add_ramping_MorLatRam13(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + T = model.instance.time + vars, eqs = model.vars, model.eqs + for t = 1:T + for g in model.instance.units + gi = g.name + known_initial_conditions = (g.initial_status != nothing && g.initial_power != nothing) + is_initially_on = known_initial_conditions && (g.initial_status > 0) + time_invariant = (t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true + RU = g.ramp_up_limit + RD = g.ramp_down_limit + + # Ramp up limit + if t == 1 + # Ignore ramping limits in first period if initial conditions are unknown + if known_initial_conditions && is_initially_on + eqs.ramp_up[gi, t] = + @constraint(mip, + g.min_power[t] + vars.prod_above[gi, t] + + (RESERVES_WHEN_RAMP_UP ? vars.reserve[gi, t] : 0.) + <= g.initial_power + RU) + end + else + # amk: without accounting for time-varying min power terms, + # we might get an infeasible schedule, e.g. if min_power[t-1] = 0, min_power[t] = 10 + # and ramp_up_limit = 5, the constraint (p'(t) + r(t) <= p'(t-1) + RU) + # would be satisfied with p'(t) = r(t) = p'(t-1) = 0 + # Note that if switch_on[t] = 1, then eqns (20) or (21) go into effect + if !time_invariant + # Use equation (24) instead + SU = g.startup_limit + max_prod_this_period = g.min_power[t] * vars.is_on[gi, t] + vars.prod_above[gi, t] + max_prod_this_period += ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? vars.reserve[gi, t] : 0. ) + min_prod_last_period = g.min_power[t-1] * vars.is_on[gi, t-1] + vars.prod_above[gi, t-1] + eqs.ramp_up[gi, t] = + @constraint(mip, + max_prod_this_period - min_prod_last_period + <= RU * vars.is_on[gi, t-1] + SU * vars.switch_on[gi, t]) + else + # Equation (26) in Kneuven et al. (2020) + # TODO what if RU < SU? places too stringent upper bound prod_above[gi, t] when starting up, and creates diff with (24). + eqs.ramp_up[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + + (RESERVES_WHEN_RAMP_UP ? vars.reserve[gi, t] : 0.) + - vars.prod_above[gi, t-1] + <= RU) + end # check time invariant + end # check if t = 1 or not + + # Ramp down limit + if t == 1 + if known_initial_conditions && is_initially_on + # TODO If RD < SD, or more specifically if + # min_power + RD < initial_power < SD + # then the generator should be able to shut down at time t = 1, + # but the constraint below will force the unit to produce power + eqs.ramp_down[gi, t] = + @constraint(mip, + g.initial_power + - (g.min_power[t] + vars.prod_above[gi, t]) + <= RD) + end + else + # amk: similar to ramp_up, need to account for time-dependent min_power + if !time_invariant + # Revert to (25) + SD = g.shutdown_limit + max_prod_last_period = g.min_power[t-1] * vars.is_on[gi, t-1] + vars.prod_above[gi, t-1] + max_prod_last_period += ( RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? vars.reserve[gi, t-1] : 0. ) # amk added + min_prod_this_period = g.min_power[t] * vars.is_on[gi, t] + vars.prod_above[gi, t] + eqs.ramp_down[gi, t] = + @constraint(mip, + max_prod_last_period - min_prod_this_period + <= RD * vars.is_on[gi, t] + SD * vars.switch_off[gi, t]) + else + # 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). + eqs.ramp_down[gi, t] = + @constraint(mip, + vars.prod_above[gi, t-1] + + (RESERVES_WHEN_RAMP_DOWN ? vars.reserve[gi, t-1] : 0.) # amk added + - vars.prod_above[gi, t] + <= RD) + end # check time invariant + end # check if t = 1 or not + end # loop over units + end # loop over time +end # add_ramping_MorLatRam13 + +# Default +Ramping_MorLatRam13 = UCComponent( + "Ramping_MorLatRam13", + "Ensure constraints on ramping are met."* + " Needs to be used in combination with shutdown rate constraints, e.g., (21b) in Kneuven et al. (2020)."* + " Based on Morales-España et al. (2013a)."* + " Eqns. (26)+(27) [replaced by (24)+(25) if time-varying min demand] in Kneuven et al. (2020).", + RampLimits, + [:is_on, :prod_above, :reserve], + [:ramp_up, :ramp_down], + add_ramping_MorLatRam13, + nothing +) # Ramping_MorLatRam13 + + +function add_ramping_DamKucRajAta16(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for g in model.instance.units + gi = g.name + known_initial_conditions = (g.initial_status != nothing && g.initial_power != nothing) + is_initially_on = known_initial_conditions && (g.initial_status > 0) + + # 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 + + for t = 1:T + time_invariant = (t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true + if t > 1 && !time_invariant + #warning("Ramping according to Damcı-Kurt et al. (2016) requires time-invariant minimum power. This does not hold for generator ", gi, ".\nmin_power[", t, "] = ", g.min_power[t], "; min_power[", t-1, "] = ", g.min_power[t-1], ". Reverting to equations (24) and (25) based on Arroyo and Conejo (2000).") + end + + max_prod_this_period = vars.prod_above[gi, t] + max_prod_this_period += ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? vars.reserve[gi, t] : 0. ) + min_prod_last_period = 0. + if t > 1 && time_invariant + min_prod_last_period = vars.prod_above[gi, t-1] + + # Equation (35) in Kneuven et al. (2020) + # Sparser version of (24) + eqs.str_ramp_up[gi, t] = + @constraint(mip, + max_prod_this_period - min_prod_last_period + <= (SU - g.min_power[t] - RU) * vars.switch_on[gi, t] + RU * vars.is_on[gi, t]) + elseif (t == 1 && known_initial_conditions && is_initially_on) || (t > 1 && !time_invariant) + # (Ignore ramping limits if initial conditions are unknown) + if t > 1 + min_prod_last_period = vars.prod_above[gi, t-1] + g.min_power[t-1] * vars.is_on[gi, t-1] + else + min_prod_last_period = max(g.initial_power, 0.) + end + + # Add the min prod at time t back in to max_prod_this_period to get _total_ production + # (instead of using the amount above minimum, as min prod for t < 1 is unknown) + max_prod_this_period += g.min_power[t] * vars.is_on[gi, t] + + # Modified version of equation (35) in Kneuven et al. (2020) + # Equivalent to (24) + eqs.str_ramp_up[gi, t] = + @constraint(mip, + max_prod_this_period - min_prod_last_period + <= (SU - RU) * vars.switch_on[gi, t] + RU * vars.is_on[gi, t]) + end + + max_prod_last_period = min_prod_last_period + max_prod_last_period += ( t > 1 && (RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN) ? vars.reserve[gi, t-1] : 0. ) # amk added + min_prod_this_period = vars.prod_above[gi, t] + on_last_period = 0. + if t > 1 + on_last_period = vars.is_on[gi, t-1] + elseif (known_initial_conditions && g.initial_status > 0) + on_last_period = 1. + end + + if t > 1 && time_invariant + # Equation (36) in Kneuven et al. (2020) + eqs.str_ramp_down[gi, t] = + @constraint(mip, + max_prod_last_period - min_prod_this_period + <= (SD - g.min_power[t] - RD) * vars.switch_off[gi, t] + RD * on_last_period) + elseif (t == 1 && known_initial_conditions && is_initially_on) || (t > 1 && !time_invariant) + # (Ignore ramping limits if initial conditions are unknown) + # Add back in min power + min_prod_this_period += g.min_power[t] * vars.is_on[gi, t] + + # Modified version of equation (36) in Kneuven et al. (2020) + # Equivalent to (25) + eqs.str_ramp_down[gi, t] = + @constraint(mip, + max_prod_last_period - min_prod_this_period + <= (SD - RD) * vars.switch_off[gi, t] + RD * on_last_period) + end + end # loop over time + end # loop over units +end # add_ramping_DamKucRajAta16 + +# Strengthened ramping +Ramping_DamKucRajAta16 = UCComponent( + "Ramping_DamKucRajAta16", + "Ensure constraints on ramping are met."* + " Based on Damcı-Kurt et al. (2016)."* + " Eqns. (35), (36) in Kneuven et al. (2020).", + RampLimits, + [:prod_above, :reserve, :is_on, :switch_on, :switch_off], + [:str_ramp_up, :str_ramp_down], + add_ramping_DamKucRajAta16, + nothing +) # Ramping_DamKucRajAta16 + + +function add_ramping_PanGua16(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for g in model.instance.units + gi = g.name + + # The following are the same for generator g across all time periods + UT = g.min_uptime + SU = g.startup_limit # startup rate, i.e., max production right after startup + SD = g.shutdown_limit # shutdown rate, i.e., max production right before shutdown + RU = g.ramp_up_limit # ramp up rate + RD = g.ramp_down_limit # ramp down rate + + for t = 1:T + Pbar = g.max_power[t] + + if Pbar < 1e-7 + # Skip this time period if max power = 0 + continue + end + + #TRD = floor((Pbar - SU) / RD) # ramp down time + # 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 Kneuven et al. (2020) + # Generalization of (20) + # Necessary that if any of the vars.switch_on = 1 in the sum, + # then vars.switch_off[gi, t+1] = 0 + eqs.str_prod_limit[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + g.min_power[t] * vars.is_on[gi, t] + + vars.reserve[gi, t] + <= Pbar * vars.is_on[gi, t] + - (t < T ? (Pbar - SD) * vars.switch_off[gi, t+1] : 0.) + - sum((Pbar - (SU + i * RU)) * vars.switch_on[gi, t-i] + for i in 0:min(UT-2, TRU, t-1)) + ) + + if UT - 2 < TRU + # Equation (40) in Kneuven et al. (2020) + # Covers an additional time period of the ramp-up trajectory, compared to (38) + eqs.prod_limit_ramp_up_extra_period[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + g.min_power[t] * vars.is_on[gi, t] + + vars.reserve[gi, t] + <= Pbar * vars.is_on[gi, t] + - sum((Pbar - (SU + i * RU)) * vars.switch_on[gi, t-i] + for i in 0:min(UT-1, TRU, t-1)) + ) + end # check UT - 2 < TRU to cover an additional time period of ramp up + + # Add in shutdown trajectory if KSD >= 0 (else this is dominated by (38)) + KSD = min( TRD, UT-1, T-t-1 ) + if KSD > 0 + KSU = min( TRU, UT-2-KSD, t-1 ) + # Equation (41) in Kneuven et al. (2020) + eqs.prod_limit_shutdown_trajectory[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + g.min_power[t] * vars.is_on[gi, t] + + (RESERVES_WHEN_SHUT_DOWN ? vars.reserve[gi, t] : 0. ) # amk added + <= Pbar * vars.is_on[gi, t] + - sum((Pbar - (SD + i * RD)) * vars.switch_off[gi, t+1+i] + for i in 0:KSD) + - sum((Pbar - (SU + i * RU)) * vars.switch_on[gi, t-i] + for i in 0:KSU) + - ((KSU >= TRU || KSU > t-2) ? 0. : + max( 0, (SU + (KSU + 1) * RU) - (SD + TRD * RD) ) * vars.switch_on[gi, t - (KSU+1)] ) + ) + end # check KSD > 0 + end # check UT > 1 + end # loop over time + end # loop over units +end # add_ramping_PanGua16 + +# Stronger variable upper bounds with ramp limits +Ramping_PanGua16 = UCComponent( + "Ramping_PanGua16", + "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 Kneuven et al. (2020)."* + " Eqns. (38), (40), (41) in Kneuven et al. (2020).", + RampLimits, + [:prod_above, :reserve, :is_on, :switch_on, :switch_off], + [:str_prod_limit, :prod_limit_ramp_up_extra_period, :prod_limit_shutdown_trajectory], + add_ramping_PanGua16, + nothing +) # Ramping_PanGua16 + + +function add_ramping_OstAnjVan12(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for g in model.instance.units + gi = g.name + + # The following are the same for generator g across all time periods + Pbar = g.max_power[t] + + if Pbar < 1e-7 + # Skip this time period if max power = 0 + continue + end + + 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 + + TRU = floor((Pbar - SU)/RU) + #TRD = floor((Pbar - SU)/RD) + # TODO check amk changed TRD wrt Kneuven et al. + TRD = ceil((Pbar - SD) / RD) # ramp down time + + # TODO check initial conditions, but maybe okay as long as (35) and (36) are also used + for t = 1:T + if UT >= 1 + # Equation (37) in Kneuven et al. (2020) + KSD = min( TRD, UT-1, T-t-1 ) + eqs.str_prod_limit[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + g.min_power[t] * vars.is_on[gi, t] + + (RESERVES_WHEN_RAMP_DOWN ? vars.reserve[gi, t] : 0.) # amk added; TODO: should this be RESERVES_WHEN_RAMP_DOWN or RESERVES_WHEN_SHUT_DOWN? + <= Pbar * vars.is_on[gi, t] + - sum((Pbar - (SD + i * RD)) * vars.switch_off[gi, t+1+i] + for i in 0:KSD) + ) + end # check UT >= 1 + end # loop over time + end # loop over units +end # add_ramping_OstAnjVan12 + + +################################################## +## Startup costs + +# TODO: should description go here, or in UCComponent.description? if here, should we eliminate UCComponent.description? +function add_startup_costs_MorLatRam13(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for g in model.instance.units + gi = g.name + S = length(g.startup_categories) + if S == 0 + continue + end + + for t in 1:T + # If unit is switching on, we must choose a startup category + # Equation (55) in Kneuven et al. (2020) + eqs.startup_choose[gi, t] = + @constraint(mip, vars.switch_on[gi, t] == sum(vars.startup[gi, s, t] for s in 1:S)) + + model.exprs.startup_cost[gi, t] = AffExpr() + + for s in 1:S + # If unit has not switched off in the last `delay` time periods, startup category is forbidden. + # The last startup category is always allowed. + if s < S + range = (t - g.startup_categories[s + 1].delay + 1):(t - g.startup_categories[s].delay) + # if initial_status < 0, then this is the amount of time the generator has been off + if g.initial_status != nothing + initial_sum = (g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0) + else + initial_sum = 0 + end + # Change of index version of equation (54) in Kneuven et al. (2020): + # startup[gi,s,t] ≤ sum_{i=s.delay}^{(s+1).delay-1} switch_off[gi,t-i] + eqs.startup_restrict[gi, s, t] = + @constraint(mip, vars.startup[gi, s, t] + <= initial_sum + sum(vars.switch_off[gi, i] for i in range if i >= 1)) + end # if s < S (not the last category) + + # Objective function terms for start-up costs + # Equation (56) in Kneuven et al. (2020) + add_to_expression!(model.exprs.startup_cost[gi, t], + vars.startup[gi, s, t], + g.startup_categories[s].cost) + end # iterate over startup categories + add_to_expression!(model.obj, model.exprs.startup_cost[gi, t]) + end # iterate over time + end # iterate over units +end # add_startup_costs_MorLatRam13 + +# Morales-España, Latorre, and Ramos, 2013 +StartupCosts_MorLatRam13 = UCComponent( + "StartupCosts_MorLatRam13", + "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 Kneuven et al. (2020)."* + " Note that the last 'constraint' is actually setting the objective."* + "\n\tstartup[gi,s,t] ≤ sum_{i=s.delay}^{(s+1).delay-1} switch_off[gi,t-i]"* + "\n\tswitch_on[gi,t] = sum_{s=1}^{length(startup_categories)} startup[gi,s,t]"* + "\n\tstartup_cost[gi,t] = sum_{s=1}^{length(startup_categories)} cost_segments[s].cost * startup[gi,s,t]", + StartupCosts, + [:startup, :switch_on, :switch_off], + [:startup_choose, :startup_restrict], + add_startup_costs_MorLatRam13, + nothing +) # StartupCosts_MorLatRam13 + + +function add_startup_costs_NowRom00(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + error("Not implemented.") +end # add_startup_costs_NowRom00 + +# Nowak and Römisch, 2000 +StartupCosts_NowRom00 = UCComponent( + "StartupCosts_NowRom00", + "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).", + StartupCosts, + nothing, + nothing, + add_startup_costs_NowRom00, + nothing +) # StartupCosts_NowRom00 + + +function add_startupcosts_KneOstWat20(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + + for g in model.instance.units + gi = g.name + S = length(g.startup_categories) + if S == 0 + continue + end + + 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 != nothing && g.initial_status < 0) ? -g.initial_status : 0 + + for t in 1:T + # 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' 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[gi, t, tmp_t], 0.; force = true) + #end + #for tmp_t in end_time+1:T + # fix(vars.downtime_arc[gi, t, tmp_t], 0.; force = true) + #end + + # Equation (59) in Kneuven 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]" + eqs.startup_at_t[gi, t] = + @constraint(mip, + vars.switch_on[gi, t] + >= sum(vars.downtime_arc[gi,tmp_t,t] + for tmp_t in t-TC+1:t-DT if tmp_t >= 1) + ) + + # Equation (60) in Kneuven et al. (2020) + # "switch_off[g,t] >= x_g(t,t') for all t' \in [t+DT, t+TC-1]" + eqs.shutdown_at_t[gi, t] = + @constraint(mip, + vars.switch_off[gi, t] + >= sum(vars.downtime_arc[gi,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 Kneuven et al. (2020) + model.exprs.startup_cost[gi, t] = AffExpr() + 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.exprs.startup_cost[gi, t], + vars.switch_on[gi, t], + g.startup_categories[default_category].cost) + + for s in 1:S-1 + # Objective function terms for start-up costs + # Equation (61) in Kneuven 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.exprs.startup_cost[gi, t], + vars.downtime_arc[gi,tmp_t,t], + g.startup_categories[s].cost - g.startup_categories[S].cost) + end + end # iterate over startup categories + add_to_expression!(model.obj, model.exprs.startup_cost[gi, t]) + end # iterate over time + end # iterate over units +end # add_startup_costs_KneOstWat20 + +function add_variables_KneOstWat20() +end # add_variables_startupcosts_KneOstWat20 + +# Kneuven, Ostrowski, and Watson, 2020 +StartupCosts_KneOstWat20 = UCComponent( + "StartupCosts_KneOstWat20", + "Extended formulation using indicator variables.", + StartupCosts, + [:switch_on, :switch_off, :downtime_arc], + [:startup_at_t, :shutdown_at_t], + add_startupcosts_KneOstWat20, + add_variables_KneOstWat20, +) # StartupCosts_KneOstWat20 + + +################################################## +## Startup / shutdown limits + +function add_startstop_limits_MorLatRam13(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for t = 1:T + for g in model.instance.units + gi = g.name + + ## 2020-10-09 amk: added eqn (20) and check of g.min_uptime + if g.min_uptime > 1 && t < T + # Equation (20) in Kneuven 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 + eqs.startstop_limit[gi,t] = + @constraint(mip, + vars.prod_above[gi, t] + vars.reserve[gi, t] + <= (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t] + - max(0, g.max_power[t] - g.startup_limit) * vars.switch_on[gi, t] + - max(0, g.max_power[t] - g.shutdown_limit) * vars.switch_off[gi, t+1]) + else + ## Startup limits + # Equation (21a) in Kneuven et al. (2020) + # Proposed by Morales-España et al. (2013a) + eqs.startup_limit[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + vars.reserve[gi, t] + <= (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t] + - max(0, g.max_power[t] - g.startup_limit) * vars.switch_on[gi, t]) + + ## Shutdown limits + if t < T + # Equation (21b) in Kneuven 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(mip, + vars.prod_above[gi, t] + + (RESERVES_WHEN_SHUT_DOWN ? vars.reserve[gi, t] : 0.) # amk added + <= (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t] + - max(0, g.max_power[t] - g.shutdown_limit) * vars.switch_off[gi, t+1]) + end + end # check if g.min_uptime > 1 + end # loop over units + end # loop over time +end # add_startstop_limits_MorLatRam13 + +StartStopLimits_MorLatRam13 = UCComponent( + "StartStopLimits_MorLatRam13", + "Startup and shutdown limits from Morales-España et al. (2013a)."* + " Eqns. (20), (21a), and (21b) in Kneuven et al. (2020).", + GenerationLimits, + [:is_on, :prod_above, :reserve, :switch_on, :switch_off], + [:startstop_limit, :startup_limit, :shutdown_limit], + add_startstop_limits_MorLatRam13, + nothing +) # StartStopLimits_MorLatRam13 + + +function add_startstop_limits_GenMorRam17(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + vars, eqs = model.vars, model.eqs + T = model.instance.time + for g in model.instance.units + gi = g.name + known_initial_conditions = (g.initial_status != nothing && g.initial_power != nothing) + + if known_initial_conditions + if g.initial_power > g.shutdown_limit + #eqs.shutdown_limit[gi, 0] = @constraint(mip, vars.switch_off[gi, 1] <= 0) + fix(vars.switch_off[gi, 1], 0.; force = true) + end + end + + for t = 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 Kneuven et al. (2020) + eqs.startstop_limit[gi,t] = + @constraint(mip, + vars.prod_above[gi, t] + vars.reserve[gi, t] + <= (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t] + - max(0, g.max_power[t] - g.startup_limit) * vars.switch_on[gi, t] + - (t < T ? max(0, g.max_power[t] - g.shutdown_limit) * vars.switch_off[gi, t+1] : 0.) + ) + else + ## Startup limits + # Equation (23a) in Kneuven et al. (2020) + eqs.startup_limit[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + vars.reserve[gi, t] + <= (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t] + - max(0, g.max_power[t] - g.startup_limit) * vars.switch_on[gi, t] + - (t < T ? max(0, g.startup_limit - g.shutdown_limit) * vars.switch_off[gi, t+1] : 0.) + ) + + ## Shutdown limits + if t < T + # Equation (23b) in Kneuven et al. (2020) + eqs.shutdown_limit[gi, t] = + @constraint(mip, + vars.prod_above[gi, t] + vars.reserve[gi, t] + <= (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t] + - (t < T ? max(0, g.max_power[t] - g.shutdown_limit) * vars.switch_off[gi, t+1] : 0.) + - max(0, g.shutdown_limit - g.startup_limit) * vars.switch_on[gi, t]) + end + end # check if g.min_uptime > 1 + end # loop over time + end # loop over units +end # add_startstop_limits_GenMorRam17 + +StartStopLimits_GenMorRam17 = UCComponent( + "StartStopLimits_GenMorRam17", + "Startup and shutdown limits from Gentile et al. (2017)."* + " Eqns. (20), (23a), and (23b) in Kneuven et al. (2020).", + GenerationLimits, + [:is_on, :prod_above, :reserve, :switch_on, :switch_off], + [:startstop_limit, :startup_limit, :shutdown_limit], + add_startstop_limits_GenMorRam17, + nothing +) # StartStopLimits_GenMorRam17 + + +################################################## +## Shutdown costs + limits + +function add_shutdown_costs_default(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + T = model.instance.time + for t = 1:T + for g in model.instance.units + gi = g.name + + shutdown_cost = 0. + if shutdown_cost > 1e-7 + # Equation (62) in Kneuven et al. (2020) + add_to_expression!(model.obj, + model.vars.switch_off[gi, t], + shutdown_cost) + end + end # loop over units + end # loop over time +end # add_shutdown_costs_default + +# Shutdown cost +DefaultShutdownCosts = UCComponent( + "DefaultShutdownCosts", + "Shutdown costs, (62) in Kneuven et al. (2020).", + ShutdownCosts, + [:switch_off], + nothing, + add_shutdown_costs_default, + nothing +) # DefaultShutdownCosts + + +################################################## +## Network constraints + +function add_network_constraints_default(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + T = model.instance.time + for lm in model.instance.lines + for t in 1:T + end # loop over time + end # loop over lines +end # add_network constraints + +DefaultNetworkConstraints = UCComponent( + "DefaultNetworkConstraints", + "Constraints for a linear approximation of the transmission network.", + SystemConstraints, + nothing, + nothing, + add_network_constraints_default, + nothing +) # DefaultNetworkConstraints + + +################################################## +## Other contraints + +function add_net_injection_default(c::UCComponent, + mip::JuMP.Model, + model::UnitCommitmentModel2) + T = model.instance.time + instance, vars, eqs, exprs = model.instance, model.vars, model.eqs, model.exprs + + if false + for t in 1:T + for b in instance.buses + # Fixed load + exprs.net_injection[b.name, t] = AffExpr(-b.load[t]) + + # Load curtailment + set_upper_bound(vars.curtail[b.name,t], b.load[t]) + add_to_expression!(exprs.net_injection[b.name, t], vars.curtail[b.name, t], 1.0) + add_to_expression!(model.obj, + vars.curtail[b.name, t], + model.instance.power_balance_penalty[t]) + end # loop over buses + + for g in instance.units + # Total production from this unit + add_to_expression!(exprs.net_injection[g.bus.name, t], vars.prod_above[g.name, t], 1.0) + add_to_expression!(exprs.net_injection[g.bus.name, t], vars.is_on[g.name, t], g.min_power[t]) + end # loop over units + + for ps in instance.price_sensitive_loads + set_upper_bound(vars.loads[ps.name, t], ps.demand[t]) + add_to_expression!(exprs.net_injection[ps.bus.name, t], vars.loads[ps.name, t], -1.0) + add_to_expression!(model.obj, vars.loads[ps.name, t], -ps.revenue[t]) + end # loop over price sensitive loads + + # Finally, now that the expression is set up, loop over buses and add the net_injection constraint + for b in instance.buses + model.eqs.net_injection_def[t, b.name] = + @constraint(model.mip, + vars.net_injection[b.name, t] + == exprs.net_injection[b.name, t]) + end # loop over buses again + + # Overall net flow should be 0 across all buses + model.eqs.power_balance[t] = + @constraint(mip, + sum(vars.net_injection[b.name, t] + for b in instance.buses) + == 0) + end # loop over time + end # if 0 (commented out) + + for t in 1:T + for ps in instance.price_sensitive_loads + # We can optionally produce extra to meet these price-sensitive loads, + # but importantly these are separate from the "fixed" demand at each bus + add_to_expression!(model.obj, vars.loads[ps.name, t], -ps.revenue[t]) + set_upper_bound(vars.loads[ps.name, t], ps.demand[t]) + end # loop over price sensitive loads + + for b in instance.buses + # If we fail to meet the demand, we suffer an (expensive) penalty + add_to_expression!(model.obj, + vars.curtail[b.name, t], + instance.power_balance_penalty[t]) + + set_upper_bound(vars.curtail[b.name,t], b.load[t]) + + # At each bus, it holds that + # total fixed demand + # = [ total production ] - [ production allocated to price-sensitve loads ] + [ slack ] + model.eqs.net_injection_def[t, b.name] = + @constraint(mip, + vars.curtail[b.name,t] + + sum(g.min_power[t] * vars.is_on[g.name,t] + vars.prod_above[g.name,t] + for g in model.instance.units + if g.bus.name == b.name) + - sum(vars.loads[ps.name,t] + for ps in instance.price_sensitive_loads + if ps.bus.name == b.name) + - b.load[t] + == vars.net_injection[b.name, t]) + end # loop over buses + + # Overall net flow should be 0 across all buses + model.eqs.power_balance[t] = + @constraint(mip, + sum(vars.net_injection[b.name, t] + for b in instance.buses) + == 0) + end # loop over time +end # add_net_injection_default + +# flow balance and price-sensitive loads +DefaultNetInjection = UCComponent( + "DefaultNetInjection", + "Ensure demand is met.", + SystemConstraints, + [:curtail, :loads, :net_injection], + [:net_injection_def, :power_balance], + add_net_injection_default, + nothing +) # DefaultNetInjection + + +################################################## +## Helpers +""" +Get union of all required variables in the components in `comps`. +""" +function get_required_variables(comps::Array{UCComponent}) :: Array{Symbol} + vars = Array{Symbol,1}() + for c in comps + if isnothing(c.vars) + continue + end + append!(vars, c.vars) + end # iterate over components + return unique(vars) +end # get_required_variables + + +""" +Get union of all required constraints in the components in `comps`. +""" +function get_required_constraints(comps::Array{UCComponent}) :: Array{Symbol} + constrs = Array{Symbol,1}() + for c in comps + if isnothing(c.constrs) + continue + end + append!(constrs, c.constrs) + end # iterate over components + return unique(constrs) +end # get_required_constraints + + +################################################## +## FORMULATIONS + +#### Default formulation #### +""" +With references to Kneuven et al. (2020) where appropriate. + +=== +Variables: + * u_g = :is_on + * v_g = :switch_on + * w_g = :switch_off + * p'_g = :prod_above + * p_g^l = :segprod + * r_g = :reserve + * s_R = :reserve_shortfall + * s_n = :curtail (any of the fixed demand not met) + * 𝛿_g^s = :startup --> different from tight formulation + * :loads --> production to meet price-sensitive demand + * :net_injection --> needed in enforce_transmission, created in DefaultNetInjection + (for a particular bus, total production - prod alloc to price-sensitive demand - total fixed demand) + * :overflow --> needed in enforce_transmission, created in DefaultRequiredComponent + * :flow --> created in enforce_transmission + +=== +Constraints: + * Uptime/Downtime: + * (2) = :binary_link, + * (3a)+(4) = :min_uptime, + * (3b)+(5) = :min_downtime, + * (?) = :switch_on_off + * Generation limits: + * (18) = :prod_limit, + * (20) = :startstop_limit, + * (21a) = :startup_limit, + * (21b) = :shutdown_limit + * Ramp limits: + * (26) = :ramp_up, + * (27) = :ramp_down + * Piecewise production: + * (42) = :segprod_limit, + * (43) = :prod_above_def, + * (44) = add to :obj + * Startup cost: + * (54) = :startup_restrict, + * (55) = :startup_choose, + * (56) = add to :obj + * Shutdown cost: + * (62) = add to :obj + * System constraints: + * (68) = :min_reserve + * kind of (65) = :power_balance +""" +DefaultFormulation = Vector{UCComponent}( + [ + DefaultRequiredConstraints, # (17) + DefaultSystemConstraints, # currently empty + DefaultGenerationLimits, # (18) + PiecewiseProduction_Garver62, # (42), (43), (44) + DefaultUpDownTime, # (2), (3), (4), (5), (?) = :switch_on_off + DefaultReserves, # (68) + Ramping_MorLatRam13, # (26), (27) + StartStopLimits_MorLatRam13, # (20), (21a), (21b) + StartupCosts_MorLatRam13, # (54), (55), (56) + DefaultShutdownCosts, # (62) + DefaultNetInjection, # kind of (65) = :power_balance + ] +) # DefaultFormulation + +""" +Same as DefaultFormulation but with (45), a weaker version of constraints (42). +""" +SparseDefaultFormulation = Vector{UCComponent}( + [ + DefaultRequiredConstraints, # (17) + DefaultSystemConstraints, # currently empty + DefaultGenerationLimits, # (18) + PiecewiseProduction_CarArr06, # (45), (43), (44) + DefaultUpDownTime, # (2), (3), (4), (5), (?) = :switch_on_off + DefaultReserves, # (68) + Ramping_MorLatRam13, # (26), (27) + StartStopLimits_MorLatRam13, # (20), (21a), (21b) + StartupCosts_MorLatRam13, # (54), (55), (56) + DefaultShutdownCosts, # (62) + DefaultNetInjection, # kind of (65) = :power_balance + ] +) # SparseDefaultFormulation + +#### Tight formulation #### +""" +From Kneuven et al. (2020), Table 3 on page 14. +Eqn :switch_on_off is used here but does not appear in the paper, I believe. + +=== +Variables using Kneuven et al. (2020) and UnitCommitment.jl notation: + * u_g = :is_on + * v_g = :switch_on + * w_g = :switch_off + * p'_g = :prod_above + * p_g^l = :segprod + * r_g = :reserve (replaces bar(p)'_g = p'_g + r_g in paper) + * s_R = :reserve_shortfall + * s_n = :curtail + * x_g = :downtime_arc --> different from default formulation + * :loads --> production to meet price-sensitive demand + * :net_injection --> needed in enforce_transmission, created in DefaultNetInjection + (for a particular bus, total production - prod alloc to price-sensitive demand - total fixed demand) + * :overflow --> needed in enforce_transmission, created in DefaultRequiredComponent + * :flow --> created in enforce_transmission + * p_{W,n} --> renewables, accounted for implicitly in :prod_above + * Θ_n --> *missing* + * f_k --> *missing* + * s_n^+ --> *missing* + * s_n^- --> *missing* + +=== +Objective: (69) = :obj + +Constraints: + * Uptime/Downtime: + * (2) = :binary_link, + * (3a)+(4) = :min_uptime, + * (3b)+(5) = :min_downtime, + * (?) = :switch_on_off + * Generation limits: + * (17) <=> nonneg on r_g = :reserve, + * (23a) = :startup_limit, + * (23b) = :shutdown_limit, + * (38) = :str_prod_limit, + * (40) if T_s^{RU} > UT_g - 2 = :prod_limit_ramp_up_extra_period, + * (41) = :prod_limit_shutdown_trajectory, + * ¿(20) = :startstop_limit? + * Ramp limits: + * (35) = :str_ramp_up, + * (36) = :str_ramp_down + * Piecewise production: + * (43) = :prod_above_def, + * (44) = add to :obj, + * (46) if UT_g > 1 = :segprod_limit, + * else (48a) and (48b) = :segprod_limita, :segprod_limitb + * Startup cost: + * (59) = :startup_at_t, + * (60) = :shutdown_at_t, + * (61) = add to :obj + * Shudown cost: (62) = add to :obj + * Network constraints: + * (64) = network constraints, handled as lazy constraints by find_violations and enforce_transmission, + * (65) = :power_balance and :net_injection_def, + * (66) = positive and negative part of slack, ** not modeled ** + * (67) = replaced by (68) = :min_reserve +""" +TightFormulation = Vector{UCComponent}( + [ + DefaultRequiredConstraints, # (17) + DefaultSystemConstraints, # currently empty + DefaultGenerationLimits, # (18) + PiecewiseProduction_KneOstWat18, # (43), (44), (46), (48) + DefaultUpDownTime, # (2), (3), (4), (5), (?) = :switch_on_off + DefaultReserves, # (68) + Ramping_DamKucRajAta16, # (35), (36) + Ramping_PanGua16, # (38), (40) + StartStopLimits_GenMorRam17, # (20), (23a), (23b) + StartupCosts_KneOstWat20, # (59), (60), (61) + DefaultShutdownCosts, # (62) + DefaultNetInjection, # take place of (65) + ] +) # TightFormulation diff --git a/src/instance.jl b/src/instance.jl index a97e9b3..3674b4c 100644 --- a/src/instance.jl +++ b/src/instance.jl @@ -8,8 +8,11 @@ using DataStructures import Base: getindex, time import GZip +abstract type UCElement end -mutable struct Bus +abstract type Time <: UCElement end + +mutable struct Bus <: UCElement name::String offset::Int load::Array{Float64} @@ -18,19 +21,19 @@ mutable struct Bus end -mutable struct CostSegment +mutable struct CostSegment <: UCElement mw::Array{Float64} cost::Array{Float64} end -mutable struct StartupCategory +mutable struct StartupCategory <: UCElement delay::Int cost::Float64 end -mutable struct Unit +mutable struct Unit <: UCElement name::String bus::Bus max_power::Array{Float64} @@ -48,10 +51,10 @@ mutable struct Unit initial_power::Union{Float64,Nothing} provides_spinning_reserves::Array{Bool} startup_categories::Array{StartupCategory} -end +end # Unit -mutable struct TransmissionLine +mutable struct TransmissionLine <: UCElement name::String offset::Int source::Bus @@ -64,19 +67,19 @@ mutable struct TransmissionLine end -mutable struct Reserves +mutable struct Reserves <: UCElement spinning::Array{Float64} end -mutable struct Contingency +mutable struct Contingency <: UCElement name::String lines::Array{TransmissionLine} units::Array{Unit} end -mutable struct PriceSensitiveLoad +mutable struct PriceSensitiveLoad <: UCElement name::String bus::Bus demand::Array{Float64} @@ -87,6 +90,8 @@ end mutable struct UnitCommitmentInstance time::Int power_balance_penalty::Array{Float64} + "Penalty for failing to meet reserve requirement." + shortfall_penalty::Array{Float64} units::Array{Unit} buses::Array{Bus} lines::Array{TransmissionLine} @@ -151,6 +156,8 @@ function from_json(json; fix=true) # Read parameters power_balance_penalty = timeseries(json["Parameters"]["Power balance penalty (\$/MW)"], default=[1000.0 for t in 1:T]) + shortfall_penalty = timeseries(json["Parameters"]["Reserve shortfall penalty (\$/MW)"], + default=[0. for t in 1:T]) # Read buses for (bus_name, dict) in json["Buses"] @@ -286,6 +293,7 @@ function from_json(json; fix=true) instance = UnitCommitmentInstance(T, power_balance_penalty, + shortfall_penalty, units, buses, lines, diff --git a/src/model2.jl b/src/model2.jl new file mode 100644 index 0000000..3cebbc2 --- /dev/null +++ b/src/model2.jl @@ -0,0 +1,475 @@ +# 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. +# Writen by Alinson S. Xavier + +using JuMP, MathOptInterface, DataStructures +import JuMP: value, fix, set_name + +# Extend some JuMP functions so that decision variables can be safely replaced by +# (constant) floating point numbers. +function value(x::Float64) + x +end + +function fix(x::Float64, v::Float64; force) + abs(x - v) < 1e-6 || error("Value mismatch: $x != $v") +end + +function set_name(x::Float64, n::String) + # nop +end + + +""" +Create a JuMP model using the variables and constraints defined by +the collection of `UCComponent`s in `formulation`. + +Parameters +=== + * `isf`: injection shift factors + * `lodf`: line outage distribution factors +""" +function build_model(; + filename::Union{String, Nothing}=nothing, + instance::Union{UnitCommitmentInstance, Nothing}=nothing, + isf::Union{Array{Float64,2}, Nothing}=nothing, + lodf::Union{Array{Float64,2}, Nothing}=nothing, + isf_cutoff::Float64=0.005, + lodf_cutoff::Float64=0.001, + optimizer=nothing, + model=nothing, + variable_names::Bool=false, + formulation::Vector{UCComponent} = UnitCommitment.DefaultFormulation, + ) :: UnitCommitmentModel2 + + if (filename == nothing) && (instance == nothing) + error("Either filename or instance must be specified") + end + + if filename != nothing + @info "Reading: $(filename)" + time_read = @elapsed begin + instance = UnitCommitment.read(filename) + end + @info @sprintf("Read problem in %.2f seconds", time_read) + end + + if length(instance.buses) == 1 + isf = zeros(0, 0) + lodf = zeros(0, 0) + else + if isf == nothing + @info "Computing injection shift factors..." + time_isf = @elapsed begin + isf = UnitCommitment.injection_shift_factors(lines=instance.lines, + buses=instance.buses) + end + @info @sprintf("Computed ISF in %.2f seconds", time_isf) + + @info "Computing line outage factors..." + time_lodf = @elapsed begin + lodf = UnitCommitment.line_outage_factors(lines=instance.lines, + buses=instance.buses, + isf=isf) + end + @info @sprintf("Computed LODF in %.2f seconds", time_lodf) + + @info @sprintf("Applying PTDF and LODF cutoffs (%.5f, %.5f)", isf_cutoff, lodf_cutoff) + isf[abs.(isf) .< isf_cutoff] .= 0 + lodf[abs.(lodf) .< lodf_cutoff] .= 0 + end + end + + @info "Building model..." + time_model = @elapsed begin + if model == nothing + if optimizer == nothing + mip = Model() + else + mip = Model(optimizer) + end + else + mip = model + end + @info "About to build model" + model = UnitCommitmentModel2(mip, # JuMP.Model + DotDict(), # vars + DotDict(), # eqs + DotDict(), # exprs + instance, # UnitCommitmentInstance + isf, # injection shift factors + lodf, # line outage distribution factors + AffExpr(), # obj + formulation, # formulation + ) + + # Prepare variables + for var in get_required_variables(formulation) + add_variable(mip, model, instance, UnitCommitment.var_list[var]) + end # prepare variables + + # Prepare constraints + for constr in get_required_constraints(formulation) + add_constraint(mip, model, instance, constr) + end # prepare constraints + + # Prepare expressions (in this case, affine expressions that are later used as part of constraints or objective) + # * :startup_cost => contribution to objective of startup costs + for field in [:startup_cost] #[:net_injection] + setproperty!(model.exprs, field, OrderedDict()) + end + + # Add components to mip + for c in formulation + c.add_component(c, mip, model) + end + + # Add objective function + build_obj_function!(model) + end # end timing of building model + @info @sprintf("Built model in %.2f seconds", time_model) + + if variable_names + set_variable_names!(model) + end + + return model +end # build_model + + +""" +Add a particular variable to `model.vars`. +""" +function add_variable(mip::JuMP.Model, + model::UnitCommitmentModel2, + instance::UnitCommitmentInstance, + var::UCVariable) + setproperty!(model.vars, var.name, OrderedDict()) + x = getproperty(model.vars, var.name) + if !isnothing(var.add_variable) + var.add_variable(var, x, mip, instance) + return + end + + # The following is a bit complex-looking, but the idea is ultimately straightforward + # We want to loop over the possible index values for var, + # for every dimension of var (e.g., looping over units and time) + # The OrderedDict `ind_to_field` maps a UCElement to the corresponding field name within a UnitCommitmentInstance + # NB: this can be an array of field names, such as [:x, :y], which means we want to access instance.x.y + # Furthermore, `var` has an array `indices` of UCElement values, describing which index loops over + # So all we want is to extract the _length_ of the corresponding field of `instance` + # We create a Tuple so we can feed it to CartesianIndices + fields = UnitCommitment.ind_to_field(var.indices) + num_indices = UnitCommitment.num_indices(fields) + + # There is some really complicated logic below that one day needs to be improved + # (we need to handle nested indices, and this is one way that hopefully works, but it is definitely not intuitive) + loop_primitive = UnitCommitment.loop_over_indices(UnitCommitment.get_indices_tuple(instance, fields)) + indices = UnitCommitment.get_indices(loop_primitive) # returns an array of tuples? or a unit range maybe. + + for ind in indices + # For each of the indices, check if the field corresponding to that index has a name + # Then we will index the variable by that name instead of the integer + curr_tuple = Tuple(ind) + new_tuple = () + for i in 1:num_indices + curr_field = UnitCommitment.get_nested_field(instance, fields, i, curr_tuple) + if :name in propertynames(curr_field) + new_tuple = (new_tuple..., curr_field.name) + else + new_tuple = (new_tuple..., curr_tuple[i]) + end + end + name = string(var.name, "[") + for (i,val) in enumerate(new_tuple) + name = string(name, val, i < num_indices ? "," : "") + end + name = string(name, "]") + if num_indices == 1 + new_tuple = new_tuple[1] + end + x[new_tuple] = @variable(mip, + lower_bound=var.lb, + upper_bound=var.ub, + integer=var.integer, + base_name=name) + end + ### DEBUG + #if var.name == :reserve_shortfall + # @show var.name, num_indices, loop_primitive, indices, x + # #@show JuMP.all_variables(mip) + #end + ### DEBUG +end # add_variable + + +""" +Add constraint to `model.eqs` (set of affine expressions represent left-hand side of constraints). +""" +function add_constraint(mip::JuMP.Model, + model::UnitCommitmentModel2, + instance::UnitCommitmentInstance, + constr::Symbol) + setproperty!(model.eqs, constr, OrderedDict()) +end # add_constraint + + +""" +Components of the objective include, summed over time: + * production cost above minimum + * minimum production cost if generator is on + * startup cost + * shutdown cost + * cost of not meeting shortfall + * penalty for not meeting or exceeding load (using curtai variable) + * shutdown cost +""" +function build_obj_function!(model::UnitCommitmentModel2) + @objective(model.mip, Min, model.obj) +end # build_obj_function + + +function enforce_transmission(; + model::UnitCommitmentModel2, + violation::Violation, + isf::Array{Float64,2}, + lodf::Array{Float64,2})::Nothing + + instance, mip, vars = model.instance, model.mip, model.vars + limit::Float64 = 0.0 + + if violation.outage_line == nothing + limit = violation.monitored_line.normal_flow_limit[violation.time] + @info @sprintf(" %8.3f MW overflow in %-5s time %3d (pre-contingency)", + violation.amount, + violation.monitored_line.name, + violation.time) + else + limit = violation.monitored_line.emergency_flow_limit[violation.time] + @info @sprintf(" %8.3f MW overflow in %-5s time %3d (outage: line %s)", + violation.amount, + violation.monitored_line.name, + violation.time, + violation.outage_line.name) + end + + fm = violation.monitored_line.name + t = violation.time + flow = @variable(mip, base_name="flow[$fm,$t]") + + # |flow| <= limit + overflow + overflow = vars.overflow[violation.monitored_line.name, violation.time] + @constraint(mip, flow <= limit + overflow) + @constraint(mip, -flow <= limit + overflow) + + if violation.outage_line == nothing + @constraint(mip, flow == sum(vars.net_injection[b.name, violation.time] * + isf[violation.monitored_line.offset, b.offset] + for b in instance.buses + if b.offset > 0)) + else + @constraint(mip, flow == sum(vars.net_injection[b.name, violation.time] * ( + isf[violation.monitored_line.offset, b.offset] + ( + lodf[violation.monitored_line.offset, violation.outage_line.offset] * + isf[violation.outage_line.offset, b.offset] + ) + ) + for b in instance.buses + if b.offset > 0)) + end + nothing +end # enforce_transmission + + +function set_variable_names!(model::UnitCommitmentModel2) + @info "Setting variable and constraint names..." + time_varnames = @elapsed begin + #set_jump_names!(model.vars) # amk: already set + set_jump_names!(model.eqs) + end + @info @sprintf("Set names in %.2f seconds", time_varnames) +end # set_variable_names + + +function set_jump_names!(dict) + for name in keys(dict) + for idx in keys(dict[name]) + idx_str = isa(idx, Tuple) ? join(map(string, idx), ",") : idx + set_name(dict[name][idx], "$name[$idx_str]") + end + end +end # set_jump_names + + +function get_solution(model::UnitCommitmentModel2) + instance, T = model.instance, model.instance.time + function timeseries(vars, collection) + return OrderedDict(b.name => [round(value(vars[b.name, t]), digits=5) for t in 1:T] + for b in collection) + end + function production_cost(g) + return [value(model.vars.is_on[g.name, t]) * g.min_power_cost[t] + + sum(Float64[value(model.vars.segprod[g.name, k, t]) * g.cost_segments[k].cost[t] + for k in 1:length(g.cost_segments)]) + for t in 1:T] + end + function production(g) + return [value(model.vars.is_on[g.name, t]) * g.min_power[t] + + sum(Float64[value(model.vars.segprod[g.name, k, t]) + for k in 1:length(g.cost_segments)]) + for t in 1:T] + end + function startup_cost(g) + #S = length(g.startup_categories) + #return [sum(g.startup_categories[s].cost * value(model.vars.startup[g.name, s, t]) + # for s in 1:S) + # for t in 1:T] + return [ value.(model.exprs.startup_cost[g.name, t]) for t in 1:T ] + end + sol = OrderedDict() + sol["Production (MW)"] = OrderedDict(g.name => production(g) for g in instance.units) + sol["Production cost (\$)"] = OrderedDict(g.name => production_cost(g) for g in instance.units) + sol["Startup cost (\$)"] = OrderedDict(g.name => startup_cost(g) for g in instance.units) + sol["Is on"] = timeseries(model.vars.is_on, instance.units) + sol["Switch on"] = timeseries(model.vars.switch_on, instance.units) + sol["Switch off"] = timeseries(model.vars.switch_off, instance.units) + sol["Reserve (MW)"] = timeseries(model.vars.reserve, instance.units) + sol["Net injection (MW)"] = timeseries(model.vars.net_injection, instance.buses) + sol["Load curtail (MW)"] = timeseries(model.vars.curtail, instance.buses) + if !isempty(instance.lines) + sol["Line overflow (MW)"] = timeseries(model.vars.overflow, instance.lines) + end + if !isempty(instance.price_sensitive_loads) + sol["Price-sensitive loads (MW)"] = timeseries(model.vars.loads, instance.price_sensitive_loads) + end + return sol +end # get_solution + + +function fix!(model::UnitCommitmentModel2, solution)::Nothing + vars, instance, T = model.vars, model.instance, model.instance.time + for g in instance.units + for t in 1:T + is_on = round(solution["Is on"][g.name][t]) + production = round(solution["Production (MW)"][g.name][t], digits=5) + reserve = round(solution["Reserve (MW)"][g.name][t], digits=5) + JuMP.fix(vars.is_on[g.name, t], is_on, force=true) + JuMP.fix(vars.prod_above[g.name, t], production - is_on * g.min_power[t], force=true) + JuMP.fix(vars.reserve[g.name, t], reserve, force=true) + end + end +end # fix! + + +function set_warm_start!(model::UnitCommitmentModel2, solution)::Nothing + vars, instance, T = model.vars, model.instance, model.instance.time + for g in instance.units + for t in 1:T + JuMP.set_start_value(vars.is_on[g.name, t], solution["Is on"][g.name][t]) + JuMP.set_start_value(vars.switch_on[g.name, t], solution["Switch on"][g.name][t]) + JuMP.set_start_value(vars.switch_off[g.name, t], solution["Switch off"][g.name][t]) + end + end +end # set_warm_start + + +function optimize!(model::UnitCommitmentModel2; + time_limit=3600, + gap_limit=1e-4, + two_phase_gap=true, + )::Nothing + + function set_gap(gap) + try + JuMP.set_optimizer_attribute(model.mip, "MIPGap", gap) + @info @sprintf("MIP gap tolerance set to %f", gap) + catch + @warn "Could not change MIP gap tolerance" + end + end + + instance = model.instance + initial_time = time() + + large_gap = false + has_transmission = (length(model.isf) > 0) + + if has_transmission && two_phase_gap + set_gap(1e-2) + large_gap = true + else + set_gap(gap_limit) + end + + while true + time_elapsed = time() - initial_time + time_remaining = time_limit - time_elapsed + if time_remaining < 0 + @info "Time limit exceeded" + break + end + + @info @sprintf("Setting MILP time limit to %.2f seconds", time_remaining) + JuMP.set_time_limit_sec(model.mip, time_remaining) + + @info "Solving MILP..." + JuMP.optimize!(model.mip) + + has_transmission || break + + violations = find_violations(model) + if isempty(violations) + @info "No violations found" + if large_gap + large_gap = false + set_gap(gap_limit) + else + break + end + else + enforce_transmission(model, violations) + end + end + + nothing +end # optimize! + + +""" +Identify which transmission lines are violated. +See find_violations description from screening.jl. +""" +function find_violations(model::UnitCommitmentModel2) + instance, vars = model.instance, model.vars + length(instance.buses) > 1 || return [] + violations = [] + @info "Verifying transmission limits..." + time_screening = @elapsed begin + non_slack_buses = [b for b in instance.buses if b.offset > 0] + net_injections = [value(vars.net_injection[b.name, t]) + for b in non_slack_buses, t in 1:instance.time] + overflow = [value(vars.overflow[lm.name, t]) + for lm in instance.lines, t in 1:instance.time] + violations = UnitCommitment.find_violations(instance=instance, + net_injections=net_injections, + overflow=overflow, + isf=model.isf, + lodf=model.lodf) + end + @info @sprintf("Verified transmission limits in %.2f seconds", time_screening) + return violations +end # find_violations + + +function enforce_transmission(model::UnitCommitmentModel2, violations::Array{Violation, 1}) + for v in violations + enforce_transmission(model=model, + violation=v, + isf=model.isf, + lodf=model.lodf) + end +end # enforce_transmission + + +export UnitCommitmentModel2, build_model, get_solution, optimize! diff --git a/src/validate.jl b/src/validate.jl index 4f0b10c..a642b58 100644 --- a/src/validate.jl +++ b/src/validate.jl @@ -115,6 +115,7 @@ function validate_units(instance, solution; tol=0.01) actual_production_cost = solution["Production cost (\$)"][unit.name] actual_startup_cost = solution["Startup cost (\$)"][unit.name] is_on = bin(solution["Is on"][unit.name]) + switch_off = bin(solution["Switch off"][unit.name]) # some formulations may not use this for t in 1:instance.time # Auxiliary variables @@ -127,7 +128,8 @@ function validate_units(instance, solution; tol=0.01) is_starting_up = !is_on[t-1] && is_on[t] is_shutting_down = is_on[t-1] && !is_on[t] ramp_up = max(0, production[t] + reserve[t] - production[t-1]) - ramp_down = max(0, production[t-1] - production[t]) + #ramp_down = max(0, production[t-1] - production[t]) + ramp_down = max(0, production[t-1] + reserve[t-1] - production[t]) end # Compute production costs @@ -193,8 +195,12 @@ function validate_units(instance, solution; tol=0.01) # Shutdown limit if is_shutting_down && (ramp_down > unit.shutdown_limit + tol) - @error @sprintf("Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)", - unit.name, t, ramp_down, unit.shutdown_limit) + @error @sprintf("Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)\n\tproduction[t-1] = %.2f\n\treserve[t-1] = %.2f\n\tproduction[t] = %.2f\n\treserve[t] = %.2f\n\tis_on[t-1] = %d\n\tis_on[t] = %d", + unit.name, t, ramp_down, unit.shutdown_limit, + (t == 1 ? unit.initial_power : production[t-1]), production[t], + (t == 1 ? 0. : reserve[t-1]), reserve[t], + (t == 1 ? unit.initial_status != nothing && unit.initial_status > 0 : is_on[t-1]), is_on[t] + ) err_count += 1 end @@ -207,8 +213,12 @@ function validate_units(instance, solution; tol=0.01) # Ramp-down limit if !is_starting_up && !is_shutting_down && (ramp_down > unit.ramp_down_limit + tol) - @error @sprintf("Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)", - unit.name, t, ramp_down, unit.ramp_down_limit) + @error @sprintf("Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)\n\tproduction[t-1] = %.2f\n\treserve[t-1] = %.2f\n\tproduction[t] = %.2f\n\treserve[t] = %.2f\n\tis_on[t-1] = %d\n\tis_on[t] = %d", + unit.name, t, ramp_down, unit.ramp_down_limit, + (t == 1 ? unit.initial_power : production[t-1]), production[t], + (t == 1 ? 0. : reserve[t-1]), reserve[t], + (t == 1 ? unit.initial_status != nothing && unit.initial_status > 0 : is_on[t-1]), is_on[t] + ) err_count += 1 end @@ -218,13 +228,16 @@ function validate_units(instance, solution; tol=0.01) # Calculate how much time the unit has been offline time_down = 0 for k in 1:(t-1) - if !is_on[t - k] + if !is_on[t - k] time_down += 1 else break end end - if t == time_down + 1 + if t == time_down + 1 && !switch_off[1] + # If unit has always been off, then the correct startup cost depends on how long was it off before t = 1 + # Absent known initial conditions, we assume it was off for the minimum downtime + # TODO: verify the formulations are making the same assumption... initial_down = unit.min_downtime if unit.initial_status < 0 initial_down = -unit.initial_status diff --git a/src/variables.jl b/src/variables.jl new file mode 100644 index 0000000..a8e163d --- /dev/null +++ b/src/variables.jl @@ -0,0 +1,471 @@ +using DataStructures # for OrderedDict +using JuMP + +################################################## +# Variables +#mutable struct UCVariable +# "Name of the variable." +# name::Symbol +# "What does the variable represent?" +# description::String +# "Global lower bound for the variable (may be adjusted later)." +# lb::Float64 +# "Global upper bound for the variable (may be adjusted later)." +# ub::Float64 +# "Is the variable integer-restricted?" +# integer::Bool +# "What are we indexing over?"* +# " Recursive structure, e.g., [X,Y] means Y is a field in X,"* +# " and [X,[Y1,Z],Y2] means Y1 and Y2 are fields in X and Z is a field in Y1.\n"* +# " [ X, [Y,A,B], [Y,A,A], [Z,[D,E],F], T ]\n"* +# " => [x, y1, y1.a, y1.b, y2, y2.a1, y2.a2, z, z.d, z.d.e, z.f, t]." +# indices::Vector +# "Function to add the variable; if this is missing, we will attempt to add the variable automatically using the `indices`. Signature should be (variable, model.vars.familyname, mip, instance)." +# add_variable::Union{Function,Nothing} +#end # UCVariable + +# TODO Above did not work for some reason +mutable struct UCVariable + name::Symbol + description::String + lb::Float64 + ub::Float64 + integer::Bool + indices::Vector + add_variable::Union{Function,Nothing} +end + +""" +It holds that 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' if the generator starts afterwards, always has max cost. +""" +function add_downtime_arcs(var::UCVariable, + x::OrderedDict, + mip::JuMP.Model, + instance::UnitCommitmentInstance) + T = instance.time + for g in instance.units + S = length(g.startup_categories) + if S == 0 + continue + end + + DT = g.min_downtime # minimum time offline + TC = g.startup_categories[S].delay # time offline until totally cold + + for t1 = 1:T-1 + for t2 = t1+1:T + # It holds that 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' if the generator starts afterwards, always has max cost + if (t2 < t1 + DT) || (t2 >= t1 + TC) + continue + end + + name = string(var.name, "[", g.name, ",", t1, ",", t2, "]") + x[g.name, t1, t2] = @variable(mip, + lower_bound=var.lb, + upper_bound=var.ub, + integer=var.integer, + base_name=name) + end # loop over time 2 + end # loop over time 1 + end # loop over units +end # add_downtime_arcs + + +""" +If there is a penalty specified for not meeting the reserve, then we add a reserve shortfall variable. +""" +function add_reserve_shortfall(var::UCVariable, + x::OrderedDict, + mip::JuMP.Model, + instance::UnitCommitmentInstance) + T = instance.time + for t = 1:T + if instance.shortfall_penalty[t] > 1e-7 + name = string(var.name, "[", t, "]") + x[t] = @variable(mip, + lower_bound=var.lb, + upper_bound=var.ub, + integer=var.integer, + base_name=name) + end + end # loop over time +end # add_reserve_shortfall + + +""" +Variables that the model may (or may not) use. + +Note the relationship + r_g(t) = bar{p}_g(t) - p_g(t) + = bar{p}'_g(t) - p'_g(t) +""" +var_list = OrderedDict{Symbol,UCVariable}( + :prod + => UCVariable(:prod, + "[gen, t]; power from generator gen at time t; p_g(t) = p'_g(t) + g.min_power[t] * u_g(t)", + 0., Inf, false, + [Unit, Time], nothing), + :prod_above + => UCVariable(:prod_above, + "[gen, t]; production above minimum required level; p'_g(t)", + 0., Inf, false, + [Unit, Time], nothing ), + :max_power_avail + => UCVariable(:max_power_avail, + "[gen, t]; maximum power available from generator gen at time t; bar{p}_g(t) = p_g(t) + r_g(t)", + 0., Inf, false, + [Unit, Time], nothing), + :max_power_avail_above + => UCVariable(:max_power_avail_above, + "[gen, t]; maximum power available above minimum from generator gen at time t; bar{p}'_g(t)", + 0., Inf, false, + [Unit, Time], nothing), + :segprod + => UCVariable(:segprod, + "[gen, seg, t]; how much generator gen produces on segment seg in time t; p_g^l(t)", + 0., Inf, false, + [ [Unit, CostSegment], Time], nothing), + :reserve + => UCVariable(:reserve, + "[gen, t]; reserves provided by gen at t; r_g(t)", + 0., Inf, false, + [Unit, Time], nothing), + :reserve_shortfall + => UCVariable(:reserve_shortfall, + "[t]; reserve shortfall at gen at t; s_R(t)", + 0., Inf, false, + [Time], add_reserve_shortfall), + :is_on + => UCVariable(:is_on, + "[gen, t]; is gen on at t; u_g(t)", + 0., 1., true, + [Unit, Time], nothing), + :switch_on + => UCVariable(:switch_on, + "[gen, t]; indicator that gen will be turned on at t; v_g(t)", + 0., 1., true, + [Unit, Time], nothing), + :switch_off + => UCVariable(:switch_off, + "[gen, t]; indicator that gen will be turned off at t; w_g(t)", + 0., 1., true, + [Unit, Time], nothing), + :net_injection + => UCVariable(:net_injection, + "[bus.name, t]", + -1e100, Inf, false, + [Bus, Time], nothing), + :curtail + => UCVariable(:curtail, + "[bus.name, t]; upper bound is max load at the bus at time t", + 0., Inf, false, + [Bus, Time], nothing), + :flow + => UCVariable(:flow, + "[violation.monitored_line.name, t]", + -1e100, Inf, false, + [Violation, Time], nothing), + :overflow + => UCVariable(:overflow, + "[transmission_line.name, t]; how much flow above the transmission limits (in MW) is allowed", + 0., Inf, false, + [TransmissionLine, Time], nothing), + :loads + => UCVariable(:loads, + "[price_sensitive_load.name, t]; production to meet demand at a set price, if it is economically sensible, independent of the rest of the demand; upper bound is demand at this price at time t", + 0., Inf, false, + [PriceSensitiveLoad, Time], nothing), + :startup + => UCVariable(:startup, + "[gen, startup_category, t]; indicator that generator g starts up in startup_category at time t; 𝛿_g^s(t)", + 0., 1., true, + [ [Unit, StartupCategory], Time], nothing), + :downtime_arc + => UCVariable(:downtime_arc, + "[gen, t, t']; indicator for shutdown at t and starting at t'", + 0., 1., true, + [Unit, Time, Time], add_downtime_arcs), +) # var_list + +#var_symbol_list = +# [ +# :prod_above, # [gen, t], ≥ 0 +# :segprod, # [gen, t, segment], ≥ 0 +# :reserve, # [gen, t], ≥ 0 +# :is_on, # [gen, t], binary +# :switch_on, # [gen, t], binary +# :switch_off, # [gen, t], binary +# :net_injection, # [bus.name, t], urs? +# :curtail, # [bus.name, t], domain [0, b.load[t]] +# :overflow, # [transmission_line.name, t], ≥ 0 +# :loads, # [price_sensitive_load.name, t], domain [0, ps.demand[t]] +# :startup # [gen, t, startup_category], binary +# ] + + +""" +For a particular UCElement, which is the field in UnitCommitmentInstance that this corresponds to? +This is used to determine indexing and ranges, e.g., `is_on` is indexed over Unit and Time, +so the variable `is_on` will range in the first index from 1 to length(instance.units) +and on the second index from 1 to instance.time. +""" +ind_to_field_dict = OrderedDict{Type{<:UCElement},Symbol}( + Time => :time, + Bus => :buses, + Unit => :units, + TransmissionLine => :lines, + PriceSensitiveLoad => :price_sensitive_loads, + CostSegment => :cost_segments, + StartupCategory => :startup_categories, +) # ind_to_field_dict + +""" +Take indices and convert them to fields of UnitCommitmentInstance. +""" +function ind_to_field(index::Union{Vector,Type{<:UCElement}}) :: Union{Vector,Symbol} + if isa(index, Type{<:UCElement}) + return ind_to_field_dict[index] + else + return [ ind_to_field(t) for t in index ] + end +end # ind_to_field + +function num_indices(v) :: Int64 + if !isa(v, Array) + return 1 + else + return sum(num_indices(v[i]) for i in 1:length(v)) + end +end # num_indices + + +""" +Can return + * UnitRange -> iterate over this range + * Array{UnitRange} -> cross product of the ranges in the array + * Tuple(UnitRange, Array{UnitRange}) -> the array length should be the same as the range of the UnitRange +""" +function get_indices_tuple(obj::Any, fields::Union{Symbol,Vector,Nothing} = nothing) + if isa(fields, Symbol) + return get_indices_tuple(getfield(obj,fields)) + end + if fields == nothing || (isa(fields,Array) && length(fields) == 0) + if isa(obj, Array) + return UnitRange(1,length(obj)) + elseif isa(obj, Int) + return UnitRange(1,obj) + else + return UnitRange{Int64}(0:-1) + #return UnitRange(1,1) + end + end + + if isa(obj,Array) + indices = ( + UnitRange(1,length(obj)), + ([ + isa(f,Array) ? get_indices_tuple(getfield(x, f[1]), f[2:end]) : get_indices_tuple(getfield(x, f)) + for x in obj + ] for f in fields)... + ) + # more_indices = ([ + # isa(f,Array) ? get_indices_tuple(getfield(x, f[1]), f[2:end]) : get_indices_tuple(getfield(x, f)) + # for x in obj + # ] for f in fields + # ) + # indices = (UnitRange(1,length(obj)),more_indices...) + else + indices = () + for f in fields + if isa(f,Array) + indices = (indices..., get_indices_tuple(getfield(obj, f[1]), f[2:end])) + else + indices = (indices..., get_indices_tuple(obj,f)) + end + end + # indices = ( + # isa(f,Array) ? get_indices_tuple(getfield(obj, f[1]), f[2:end]) : get_indices_tuple(getfield(obj, f)) + # for f in fields + # ) + # ( + # isa(f,Array) ? get_indices_tuple(getfield(obj, f[1]), f[2:end]) : get_indices_tuple(getfield(obj, f)) + # for f in fields + # ) + # indices = (indices...,) + end # check if obj is Array or not + + return indices +end # get_indices_tuple + +function loop_over_indices(indices::Any) + loop = nothing + should_print = false + + if isa(indices, UnitRange) + loop = indices + elseif isa(indices, Array{UnitRange{Int64}}) || isa(indices, Tuple{Int, UnitRange}) + loop = Base.product(Tuple(indices)...) + elseif isa(indices, Tuple{UnitRange, Array}) + loop = () + for t in zip(indices...) + loop = (loop..., loop_over_indices(t)...) + end + elseif isa(indices,Tuple) + loop = () + for i in indices + loop = (loop..., loop_over_indices(i)) + end + loop = Base.product(loop...) + else + error("Why are we here?") + #loop = Base.product(loop_over_indices(indices)...) + end + + if should_print + for i in loop + @show i + end + end + return loop +end # loop_over_indices + + +function expand_tuple(x::Tuple) + y = () + for i in x + if isa(i, Tuple) + y = (y..., expand_tuple(i)...) + else + y = (y..., i) + end + end + return y +end # expand_tuple + + +function expand_tuple(X::Array{<:Tuple}) + return [ expand_tuple(x) for x in X ] +end # expand_tuple + + +function get_indices(x::Array) + return expand_tuple(x) +end +function get_indices(x::Base.Iterators.ProductIterator) + return get_indices(collect(x)) +end + + +""" +Access `t.f`, special terminal case of `get_nested_field`. +""" +function get_nested_field(t::Any, f::Symbol) + return getfield(t,f) +end # get_nested_field + + +""" +Access `t.f`, where `f` could be a subfield. +""" +function get_nested_field(t::Any, f::Vector{Symbol}) + if length(f) > 1 + return get_nested_field(getfield(t,f[1]), f[2:end]) + else + return getfield(t,f[1]) + end +end # get_nested_field + + +""" +Given a set of indices of UCVariable, e.g., [[X,Y],T], +and a UnitCommitmentInstance instance, +if we want to access the field corresponding to Y, +then we call get_nested_field(instance, [[X,Y],T], 2, (4,3)), +which will return instance.X[4].Y[3] if Y is a vector, +and just instance.X[4].Y otherwise. + +=== +Termination Conditions + +If `i` <= 0, then we only care about instance, and not the field. +If `field` is a Symbol and `i` >= 1, then we want to explore instance.field (index t[i] or t). +Note that if `i` >= 1, then `field` must be a symbol. +If `i` == 1, then `t` can be an Int. + +=== +Parameters + * instance::Any --> all fields will be from instance, or nested fields of fields of instance. + * field::Union{Vector,Symbol,Nothing} --> either the field we want to access, or a vector of fields, and we will want field[i]. + * i::Int --> which field to access. + * t::Tuple --> how to go through the fields of instance to get the right field, length needs to be at least `i`. +""" +function get_nested_field(instance::Any, field::Union{Vector,Symbol,Nothing}, i::Int, t::Union{Tuple, Int}) + # Check data + if isa(field, Vector) + if i >= 2 && (!isa(t,Tuple) || length(t) < i) + error("Tuple of indices to get nested field needs to be at least the length of the index we want to get.") + end + end + + if isa(field, Symbol) || i <= 0 + # i = 0 can happen in the recursive call + # What it means is that we do not want a field of the instance, but the instance itself + # TODO handle other iterable types and empty arrays + f = (isa(field, Symbol) && i >= 1) ? getfield(instance, field) : instance + if isa(f,Vector) + if length(f) == 0 + error("Trying to iterate over empty field!") + else + return isa(t,Int) ? f[t] : f[t[i]] + end + else + return f + end + end # check termination conditions (f is field or i <= 0) + + # Loop over the fields until we find where index i is located + # It may be nested inside an array, so that is why we recurse + start_ind = 0 + for f in field + curr_len = isa(f, Vector) ? length(f) : 1 + if start_ind + curr_len >= i + if isa(f, Vector) + new_field_is_iterable = isa(getfield(instance, f[1]), Vector) + if new_field_is_iterable + return get_nested_field(getfield(instance, f[1])[t[start_ind+1]], f[2:end], i - start_ind - 1, isa(t,Tuple) ? t[start_ind+2:end] : t) + else + return get_nested_field(getfield(instance, f[1]), f[2:end], i - start_ind - 1, isa(t,Tuple) ? t[start_ind+2:end] : t) + end + else + # f is hopefully a symbol... + return get_nested_field(instance, f, 1, isa(t,Tuple) ? t[start_ind+1] : t) + end + end + start_ind += curr_len + end + + return nothing +end # get_nested_field + + +#""" +#Get ranges for the indices of a UCVariable along dimension `i`, +#making sure that the right fields ranges are calculated via `get_nested_field` and `ind_to_field`. +#""" +#function get_range(arr::UCVariable, instance::UnitCommitmentInstance, i::Int) :: UnitRange +# arr = ind_to_field[var.indices[i]] +# f = get_nested_field(instance, arr) +# if isa(f, Array) +# return 1:length(f) +# elseif isa(f, Int) +# return 1:f +# else +# error("Unknown type to generate UnitRange from: ", typeof(f)) +# end +#end # get_range + +export UCVariable diff --git a/test/model_test.jl b/test/model_test.jl index 7d75d2e..f1fdddf 100644 --- a/test/model_test.jl +++ b/test/model_test.jl @@ -2,31 +2,55 @@ # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Released under the modified BSD license. See COPYING.md for more details. -using UnitCommitment, LinearAlgebra, Cbc, JuMP +using UnitCommitment, LinearAlgebra, JuMP +USE_GUROBI = (Base.find_package("Gurobi") != nothing) +USE_CBC = !USE_GUROBI +if USE_GUROBI + using Gurobi +else + using Cbc +end + +NUM_THREADS = 4 +LOG_LEVEL = 1 @testset "Model" begin @testset "Run" begin - instance = UnitCommitment.read_benchmark("test/case14") + #instance = UnitCommitment.read_benchmark("test/case14") + #instance = UnitCommitment.read_benchmark("matpower/case3375wp/2017-02-01") + instance = UnitCommitment.read_benchmark("matpower/case1888rte/2017-02-01") for line in instance.lines, t in 1:4 line.normal_flow_limit[t] = 10.0 end - optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) - model = build_model(instance=instance, - optimizer=optimizer, - variable_names=true) + #for formulation in [UnitCommitment.DefaultFormulation, UnitCommitment.TightFormulation] + for formulation in [UnitCommitment.TightFormulation] + @info string("Running test of ", formulation) + if USE_CBC + optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => LOG_LEVEL) + end + if USE_GUROBI + optimizer = optimizer_with_attributes(Gurobi.Optimizer, "Threads" => NUM_THREADS) + end + model = build_model(instance=instance, + optimizer=optimizer, + variable_names=true, + components=formulation) - JuMP.write_to_file(model.mip, "test.mps") + JuMP.write_to_file(model.mip, "test.mps") - # Optimize and retrieve solution - UnitCommitment.optimize!(model) - solution = get_solution(model) - - # Verify solution - @test UnitCommitment.validate(instance, solution) + # Optimize and retrieve solution + UnitCommitment.optimize!(model) + solution = get_solution(model) + + # Verify solution + @test UnitCommitment.validate(instance, solution) - # Reoptimize with fixed solution - UnitCommitment.fix!(model, solution) - UnitCommitment.optimize!(model) - @test UnitCommitment.validate(instance, solution) - end -end + # Reoptimize with fixed solution + UnitCommitment.fix!(model, solution) + UnitCommitment.optimize!(model) + @test UnitCommitment.validate(instance, solution) + + #@show solution + end # loop over components + end # end testset Run +end # end test