diff --git a/benchmarks/OptimizationFrameworks/Manifest.toml b/benchmarks/OptimizationFrameworks/Manifest.toml index 238571c69..662a341d0 100644 --- a/benchmarks/OptimizationFrameworks/Manifest.toml +++ b/benchmarks/OptimizationFrameworks/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.9.1" manifest_format = "2.0" -project_hash = "5356efadf00ee4c7d5ece025a6be4d6b70be771b" +project_hash = "41446dc14eec22d904f52b0fba22212fba705517" [[deps.ADNLPModels]] deps = ["ColPack", "ForwardDiff", "LinearAlgebra", "NLPModels", "Requires", "ReverseDiff", "SparseArrays"] @@ -244,6 +244,11 @@ git-tree-sha1 = "02d2316b7ffceff992f3096ae48c7829a8aa0638" uuid = "b152e2b5-7a66-4b01-a709-34e65c35f657" version = "0.1.3" +[[deps.ConcreteStructs]] +git-tree-sha1 = "f749037478283d372048690eb3b5f92a79432b34" +uuid = "2569d6c7-a4a2-43d3-a901-331e8e4be471" +version = "0.2.3" + [[deps.CondaPkg]] deps = ["JSON3", "Markdown", "MicroMamba", "Pidfile", "Pkg", "TOML"] git-tree-sha1 = "741146cf2ced5859faae76a84b541aa9af1a78bb" @@ -283,6 +288,12 @@ git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" version = "1.15.0" +[[deps.DataFrames]] +deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "REPL", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] +git-tree-sha1 = "04c738083f29f86e62c8afc341f0967d8717bdb8" +uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +version = "1.6.1" + [[deps.DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] git-tree-sha1 = "cf25ccb972fec4e4817764d01c82386ae94f77b4" @@ -558,6 +569,12 @@ git-tree-sha1 = "dc1e2eba1a98aa457b629fe1723d9078ecb74340" uuid = "2030c09a-7f63-5d83-885d-db604e0e9cc0" version = "0.7.7" +[[deps.InlineStrings]] +deps = ["Parsers"] +git-tree-sha1 = "9cc2baf75c6d09f9da536ddf58eb2f29dedaf461" +uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48" +version = "1.4.0" + [[deps.IntegerMathUtils]] git-tree-sha1 = "b8ffb903da9f7b8cf695a8bead8e01814aa24b30" uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235" @@ -579,6 +596,11 @@ git-tree-sha1 = "16c0cc91853084cb5f58a78bd209513900206ce6" uuid = "8197267c-284f-5f27-9208-e0e47529a953" version = "0.7.4" +[[deps.InvertedIndices]] +git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038" +uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f" +version = "1.3.0" + [[deps.Ipopt]] deps = ["Ipopt_jll", "LinearAlgebra", "MathOptInterface", "OpenBLAS32_jll", "PrecompileTools"] git-tree-sha1 = "fdb3430a3b7b909bcb5abf9439d61450f49e3e2c" @@ -1064,6 +1086,12 @@ git-tree-sha1 = "240d7170f5ffdb285f9427b92333c3463bf65bf6" uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad" version = "0.2.1" +[[deps.PooledArrays]] +deps = ["DataAPI", "Future"] +git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7" +uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" +version = "1.4.2" + [[deps.PositiveFactorizations]] deps = ["LinearAlgebra"] git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20" @@ -1098,6 +1126,12 @@ git-tree-sha1 = "7eb1686b4f04b82f96ed7a4ea5890a4f0c7a09f1" uuid = "21216c6a-2e73-6563-6e65-726566657250" version = "1.4.0" +[[deps.PrettyTables]] +deps = ["Crayons", "LaTeXStrings", "Markdown", "Printf", "Reexport", "StringManipulation", "Tables"] +git-tree-sha1 = "ee094908d720185ddbdc58dbe0c1cbe35453ec7a" +uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +version = "2.2.7" + [[deps.Primes]] deps = ["IntegerMathUtils"] git-tree-sha1 = "4c9f306e5d6603ae203c2000dd460d81a5251489" @@ -1251,6 +1285,12 @@ git-tree-sha1 = "30449ee12237627992a99d5e30ae63e4d78cd24a" uuid = "6c6a2e73-6563-6170-7368-637461726353" version = "1.2.0" +[[deps.SentinelArrays]] +deps = ["Dates", "Random"] +git-tree-sha1 = "04bdff0b09c65ff3e06a05e3eb7b120223da3d39" +uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c" +version = "1.4.0" + [[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -1328,6 +1368,12 @@ weakdeps = ["ChainRulesCore"] [deps.SpecialFunctions.extensions] SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" +[[deps.StableRNGs]] +deps = ["Random", "Test"] +git-tree-sha1 = "3be7d49667040add7ee151fefaf1f8c04c8c8276" +uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" +version = "1.0.0" + [[deps.Static]] deps = ["IfElse"] git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc" @@ -1400,6 +1446,12 @@ git-tree-sha1 = "f02eb61eb5c97b48c153861c72fbbfdddc607e06" uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da" version = "0.4.17" +[[deps.StringManipulation]] +deps = ["PrecompileTools"] +git-tree-sha1 = "25ad9b0ccbc6cc63dfff1eb3f4bbb4d3e0d8a491" +uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e" +version = "0.3.2" + [[deps.StructArrays]] deps = ["Adapt", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"] git-tree-sha1 = "521a0e828e98bb69042fec1809c1b5a680eb7389" diff --git a/benchmarks/OptimizationFrameworks/Project.toml b/benchmarks/OptimizationFrameworks/Project.toml index de443fbc9..279899bc0 100644 --- a/benchmarks/OptimizationFrameworks/Project.toml +++ b/benchmarks/OptimizationFrameworks/Project.toml @@ -1,6 +1,9 @@ [deps] ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" @@ -12,7 +15,9 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] diff --git a/benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m b/benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m new file mode 100644 index 000000000..b52e7555c --- /dev/null +++ b/benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m @@ -0,0 +1,103 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% %%%%% +%%%% IEEE PES Power Grid Library - Optimal Power Flow - v23.07 %%%%% +%%%% (https://github.com/power-grid-lib/pglib-opf) %%%%% +%%%% Benchmark Group - Typical Operations %%%%% +%%%% 23 - July - 2023 %%%%% +%%%% %%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% The semidefinite relaxation of the OPF problem successfully solves +% this network with a value of 60 MVA for the line-flow limit on the line from +% bus 2 to bus 3. The semidefinite relaxation fails to give a physically +% meaningful solution to this network with a value of 50 MVA for the line-flow +% limit on this line. See the following publication for further details. +% +% Lesieutre, B.C. & Molzahn, D.K. & Borden, AR. & Demarco, C.L., +% "Examining the Limits of the Application of Semidefinite Programming to Power Flow Problems", +% 49th Annual Allerton Conference on Communication, Control, and Computing (Allerton), +% September, 2011, pp. 1492-1499 +% +% opt objective value: 5812.64 $/hr +% +% Bus Voltage Generation Load Lambda($/MVA-hr) +% # Mag(pu) Ang(deg) P (MW) Q (MVAr) P (MW) Q (MVAr) P Q +% ----- ------- -------- -------- -------- -------- -------- ------- ------- +% 1 1.100 0.000* 148.07 54.70 110.00 40.00 37.575 - +% 2 0.926 7.259 170.01 -8.79 110.00 40.00 30.101 - +% 3 0.900 -17.267 0.00 -4.84 95.00 50.00 45.537 - +% -------- -------- -------- -------- +% Total: 318.07 41.06 315.00 130.00 +% +% Copyright (c) 2011 by The Institute of Electrical and Electronics Engineers (IEEE) +% Licensed under the Creative Commons Attribution 4.0 +% International license, http://creativecommons.org/licenses/by/4.0/ +% +% Contact M.E. Brennan (me.brennan@ieee.org) for inquries on further reuse of +% this dataset. +% +function mpc = pglib_opf_case3_lmbd +mpc.version = '2'; +mpc.baseMVA = 100.0; + +%% bus data +% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin +mpc.bus = [ + 1 3 110.0 40.0 0.0 0.0 1 1.00000 0.00000 240.0 1 1.10000 0.90000; + 2 2 110.0 40.0 0.0 0.0 1 1.00000 0.00000 240.0 1 1.10000 0.90000; + 3 2 95.0 50.0 0.0 0.0 1 1.00000 0.00000 240.0 1 1.10000 0.90000; +]; + +%% generator data +% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin +mpc.gen = [ + 1 1000.0 0.0 1000.0 -1000.0 1.0 100.0 1 2000.0 0.0; + 2 1000.0 0.0 1000.0 -1000.0 1.0 100.0 1 2000.0 0.0; + 3 0.0 0.0 1000.0 -1000.0 1.0 100.0 1 0.0 0.0; +]; + +%% generator cost data +% 2 startup shutdown n c(n-1) ... c0 +mpc.gencost = [ + 2 0.0 0.0 3 0.110000 5.000000 0.000000; + 2 0.0 0.0 3 0.085000 1.200000 0.000000; + 2 0.0 0.0 3 0.000000 0.000000 0.000000; +]; + +%% branch data +% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax +mpc.branch = [ + 1 3 0.065 0.62 0.45 9000.0 9000.0 9000.0 0.0 0.0 1 -30.0 30.0; + 3 2 0.025 0.75 0.7 50.0 50.0 50.0 0.0 0.0 1 -30.0 30.0; + 1 2 0.042 0.9 0.3 9000.0 9000.0 9000.0 0.0 0.0 1 -30.0 30.0; +]; + +% INFO : === Translation Options === +% INFO : Phase Angle Bound: 30.0 (deg.) +% INFO : Setting Flat Start +% INFO : +% INFO : === Generator Bounds Update Notes === +% INFO : +% INFO : === Base KV Replacement Notes === +% INFO : +% INFO : === Transformer Setting Replacement Notes === +% INFO : +% INFO : === Line Capacity Monotonicity Notes === +% INFO : Updated Thermal Rating: on line 1-3 : Rate B, Rate C - 0.0, 0.0 -> 9000.0, 9000.0 +% INFO : Updated Thermal Rating: on line 3-2 : Rate B, Rate C - 0.0, 0.0 -> 50.0, 50.0 +% INFO : Updated Thermal Rating: on line 1-2 : Rate B, Rate C - 0.0, 0.0 -> 9000.0, 9000.0 +% INFO : +% INFO : === Voltage Setpoint Replacement Notes === +% INFO : Bus 1 : V=1.0, theta=0.0 -> V=1.0, theta=0.0 +% INFO : Bus 2 : V=1.0, theta=0.0 -> V=1.0, theta=0.0 +% INFO : Bus 3 : V=1.0, theta=0.0 -> V=1.0, theta=0.0 +% INFO : +% INFO : === Generator Setpoint Replacement Notes === +% INFO : Gen at bus 1 : Pg=0.0, Qg=0.0 -> Pg=1000.0, Qg=0.0 +% INFO : Gen at bus 1 : Vg=1.0 -> Vg=1.0 +% INFO : Gen at bus 2 : Pg=0.0, Qg=0.0 -> Pg=1000.0, Qg=0.0 +% INFO : Gen at bus 2 : Vg=1.0 -> Vg=1.0 +% INFO : Gen at bus 3 : Pg=0.0, Qg=0.0 -> Pg=0.0, Qg=0.0 +% INFO : Gen at bus 3 : Vg=1.0 -> Vg=1.0 +% INFO : +% INFO : === Writing Matpower Case File Notes === diff --git a/benchmarks/OptimizationFrameworks/optimal_powerflow.jmd b/benchmarks/OptimizationFrameworks/optimal_powerflow.jmd index 102a74de3..65bf821a7 100644 --- a/benchmarks/OptimizationFrameworks/optimal_powerflow.jmd +++ b/benchmarks/OptimizationFrameworks/optimal_powerflow.jmd @@ -3,35 +3,164 @@ title: Optimal Powerflow Nonlinear Optimization Benchmark author: Chris Rackauckas --- -## Data Load +## Data Load and Setup Code + +This is generic setup code usable for all solver setups. Basically removing some unnecessary untyped dictionaries +before getting to the benchmarks. ```julia import PowerModels -time_data_start = time() +import ConcreteStructs +using BenchmarkTools + +ConcreteStructs.@concrete struct DataRepresentation + data + ref + var_lookup + var_init + var_lb + var_ub + ref_gen_idxs + lookup_pg + lookup_qg + lookup_va + lookup_vm + lookup_lij + lookup_p_lij + lookup_q_lij + cost_arrs + f_bus + t_bus + ref_bus_idxs + ref_buses_idxs + ref_bus_gens + ref_bus_arcs + ref_branch_idxs + ref_arcs_from + ref_arcs_to + p_idxmap + q_idxmap + bus_pd + bus_qd + bus_gs + bus_bs + br_g + br_b + br_tr + br_ti + br_ttm + br_g_fr + br_b_fr + br_g_to + br_b_to +end -file_name = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m" -const data = PowerModels.parse_file(file_name) -PowerModels.standardize_cost_terms!(data, order=2) -PowerModels.calc_thermal_limits!(data) -const ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] +function load_and_setup_data(file_name) + data = PowerModels.parse_file(file_name) + PowerModels.standardize_cost_terms!(data, order=2) + PowerModels.calc_thermal_limits!(data) + ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] -data_load_time = time() - time_data_start -``` + # Some data munging to type-stable forms -## Optimization.jl + var_lookup = Dict{String,Int}() -Constraint optimization implementation reference: https://github.com/SciML/Optimization.jl/blob/master/lib/OptimizationMOI/test/runtests.jl -Other AD libraries can be considered: https://docs.sciml.ai/dev/modules/Optimization/API/optimization_function/ -However ForwardDiff is the only one that is compatible with constraint functions + var_init = Float64[] + var_lb = Float64[] + var_ub = Float64[] -```julia -import Optimization -import OptimizationMOI -import ModelingToolkit -import Ipopt + var_idx = 1 + for (i,bus) in ref[:bus] + push!(var_init, 0.0) #va + push!(var_lb, -Inf) + push!(var_ub, Inf) + var_lookup["va_$(i)"] = var_idx + var_idx += 1 -function solve_opf_optimization(file_name) - time_model_start = time() + push!(var_init, 1.0) #vm + push!(var_lb, bus["vmin"]) + push!(var_ub, bus["vmax"]) + var_lookup["vm_$(i)"] = var_idx + var_idx += 1 + end + + for (i,gen) in ref[:gen] + push!(var_init, 0.0) #pg + push!(var_lb, gen["pmin"]) + push!(var_ub, gen["pmax"]) + var_lookup["pg_$(i)"] = var_idx + var_idx += 1 + + push!(var_init, 0.0) #qg + push!(var_lb, gen["qmin"]) + push!(var_ub, gen["qmax"]) + var_lookup["qg_$(i)"] = var_idx + var_idx += 1 + end + + for (l,i,j) in ref[:arcs] + branch = ref[:branch][l] + + push!(var_init, 0.0) #p + push!(var_lb, -branch["rate_a"]) + push!(var_ub, branch["rate_a"]) + var_lookup["p_$(l)_$(i)_$(j)"] = var_idx + var_idx += 1 + + push!(var_init, 0.0) #q + push!(var_lb, -branch["rate_a"]) + push!(var_ub, branch["rate_a"]) + var_lookup["q_$(l)_$(i)_$(j)"] = var_idx + var_idx += 1 + end + + @assert var_idx == length(var_init)+1 + + ref_gen_idxs = [i for i in keys(ref[:gen])] + lookup_pg = Dict{Int,Int}() + lookup_qg = Dict{Int,Int}() + lookup_va = Dict{Int,Int}() + lookup_vm = Dict{Int,Int}() + lookup_lij = Tuple{Int,Int,Int}[] + lookup_p_lij = Int[] + lookup_q_lij = Int[] + cost_arrs = Dict{Int,Vector{Float64}}() + + for (i,gen) in ref[:gen] + lookup_pg[i] = var_lookup["pg_$(i)"] + lookup_qg[i] = var_lookup["qg_$(i)"] + cost_arrs[i] = gen["cost"] + end + + for (i,bus) in ref[:bus] + lookup_va[i] = var_lookup["va_$(i)"] + lookup_vm[i] = var_lookup["vm_$(i)"] + end + + for (l,i,j) in ref[:arcs] + push!(lookup_lij, (l,i,j)) + push!(lookup_p_lij,var_lookup["p_$(l)_$(i)_$(j)"]) + push!(lookup_q_lij,var_lookup["q_$(l)_$(i)_$(j)"]) + end + + f_bus = Dict{Int,Int}() + t_bus = Dict{Int,Int}() + + for (l,branch) in ref[:branch] + f_bus[l] = branch["f_bus"] + t_bus[l] = branch["t_bus"] + end + + ref_bus_idxs = [i for i in keys(ref[:bus])] + ref_buses_idxs = [i for i in keys(ref[:ref_buses])] + ref_bus_gens = ref[:bus_gens] + ref_bus_arcs = ref[:bus_arcs] + ref_branch_idxs = [i for i in keys(ref[:branch])] + ref_arcs_from = ref[:arcs_from] + ref_arcs_to = ref[:arcs_to] + + p_idxmap = Dict(lookup_lij[i] => lookup_p_lij[i] for i in 1:length(lookup_lij)) + q_idxmap = Dict(lookup_lij[i] => lookup_q_lij[i] for i in 1:length(lookup_lij)) bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus]) bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus]) @@ -81,186 +210,247 @@ function solve_opf_optimization(file_name) br_b_to[i] = branch["b_to"] end - var_lookup = Dict{String,Int}() + DataRepresentation( + data, + ref, + var_lookup, + var_init, + var_lb, + var_ub, + ref_gen_idxs, + lookup_pg, + lookup_qg, + lookup_va, + lookup_vm, + lookup_lij, + lookup_p_lij, + lookup_q_lij, + cost_arrs, + f_bus, + t_bus, + ref_bus_idxs, + ref_buses_idxs, + ref_bus_gens, + ref_bus_arcs, + ref_branch_idxs, + ref_arcs_from, + ref_arcs_to, + p_idxmap, + q_idxmap, + bus_pd, + bus_qd, + bus_gs, + bus_bs, + br_g, + br_b, + br_tr, + br_ti, + br_ttm, + br_g_fr, + br_b_fr, + br_g_to, + br_b_to) +end - var_init = Float64[] - var_lb = Float64[] - var_ub = Float64[] +file_name = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m" +dataset = load_and_setup_data(file_name); +``` - var_idx = 1 - for (i,bus) in ref[:bus] - push!(var_init, 0.0) #va - push!(var_lb, -Inf) - push!(var_ub, Inf) - var_lookup["va_$(i)"] = var_idx - var_idx += 1 +## Test Setup - push!(var_init, 1.0) #vm - push!(var_lb, bus["vmin"]) - push!(var_ub, bus["vmax"]) - var_lookup["vm_$(i)"] = var_idx - var_idx += 1 - end +Ensure that all objectives and constraints evaluate to the same value on a random vector on the same dataset - for (i,gen) in ref[:gen] - push!(var_init, 0.0) #pg - push!(var_lb, gen["pmin"]) - push!(var_ub, gen["pmax"]) - var_lookup["pg_$(i)"] = var_idx - var_idx += 1 +```julia +test_u0 = [0.6292298794022337, 0.30740951571225206, 0.0215258802699263, 0.38457509230779996, 0.9419186480931858, 0.34961116773074874, 0.875763562401991, 0.3203478635827923, 0.6354060958226175, 0.45537545721771266, 0.3120599359696674, 0.2421238802331842, 0.886455177641366, 0.49797378087768696, 0.652913329799645, 0.03590201299300255, 0.5618806749518928, 0.8142146688533769, 0.3973557130434364, 0.27827135011662674, 0.16456134856048643, 0.7465018431665373, 0.4898329811551083, 0.6966035226583556, 0.7419662648518377, 0.8505905798503723, 0.27102126066405097, 0.1988238097281576, 0.09684601934490256, 0.49238142828542797, 0.1366594202307445, 0.6337080281764231, 0.28814906958008235, 0.5404996094640431, 0.015153517398975858, 0.6338449294034381, 0.5165464961007717, 0.572879113636733, 0.9652420600585092, 0.26535868365228543, 0.865686920119479, 0.38426996353892773, 0.007412077949221274, 0.3889835001514599] +test_obj = 7079.190664351089 +test_cons = [0.0215258802699263, -1.0701734802505833, -5.108902216849063, -3.49724505910433, -2.617834191007569, 0.5457423426033834, -0.7150251969424766, -2.473175092089014, -2.071687022809815, -1.5522321037165985, -1.0107399030803794, 3.0047739260369246, 0.2849522377447594, -2.8227966798520674, 3.2236954017592256, 1.0793383525116511, -1.633412293595111, -3.1618224299953224, -0.7775962590542184, 1.7252573527333024, -4.23535583005632, -1.7030832394691608, 1.5810450617647889, -0.33289810365419437, 0.19476447251065077, 1.0688558672739048, 1.563372246165339, 9.915310272572729, 1.4932615291788414, 2.0016715378998793, -1.4038702698147258, -0.8834081057449231, 0.21730536348839036, -7.40879932706212, -1.6000837514115611, 0.8542376821320647, 0.06615508569119477, -0.6077039991323074, 0.6138802155526912, 0.0061762164203837955, -0.3065125522705683, 0.5843454392910835, 0.7251928172073308, 1.2740182727083802, 0.11298343104675009, 0.2518186223833513, 0.4202616621130535, 0.3751697141306502, 0.4019890236200105, 0.5950107614751935, 1.0021074654956683, 0.897077248544158, 0.15136310228960612] +``` - push!(var_init, 0.0) #qg - push!(var_lb, gen["qmin"]) - push!(var_ub, gen["qmax"]) - var_lookup["qg_$(i)"] = var_idx - var_idx += 1 - end +## Setup and Validations - for (l,i,j) in ref[:arcs] - branch = ref[:branch][l] +Now is the setup code for each optimization framework, along with the validation runs on the test case. Any test which fails the validation +case, i.e. has `x_test_res[1] !≈ test_obj` or `x_test_res[2] !≈ test_cons` should be considered invalidated as this means that the model +in that modeling platform does not evaluate to give the same results - push!(var_init, 0.0) #p - push!(var_lb, -branch["rate_a"]) - push!(var_ub, branch["rate_a"]) - var_lookup["p_$(l)_$(i)_$(j)"] = var_idx - var_idx += 1 +### Optimization.jl - push!(var_init, 0.0) #q - push!(var_lb, -branch["rate_a"]) - push!(var_ub, branch["rate_a"]) - var_lookup["q_$(l)_$(i)_$(j)"] = var_idx - var_idx += 1 - end +Constraint optimization implementation reference: https://github.com/SciML/Optimization.jl/blob/master/lib/OptimizationMOI/test/runtests.jl +Other AD libraries can be considered: https://docs.sciml.ai/dev/modules/Optimization/API/optimization_function/ - @assert var_idx == length(var_init)+1 +```julia +import Optimization +import OptimizationMOI +import ModelingToolkit +import Ipopt + +function build_opf_optimization_prob(dataset; adchoice = Optimization.AutoForwardDiff()) + (;data, + ref, + var_lookup, + var_init, + var_lb, + var_ub, + ref_gen_idxs, + lookup_pg, + lookup_qg, + lookup_va, + lookup_vm, + lookup_lij, + lookup_p_lij, + lookup_q_lij, + cost_arrs, + f_bus, + t_bus, + ref_bus_idxs, + ref_buses_idxs, + ref_bus_gens, + ref_bus_arcs, + ref_branch_idxs, + ref_arcs_from, + ref_arcs_to, + p_idxmap, + q_idxmap, + bus_pd, + bus_qd, + bus_gs, + bus_bs, + br_g, + br_b, + br_tr, + br_ti, + br_ttm, + br_g_fr, + br_b_fr, + br_g_to, + br_b_to) = dataset #total_callback_time = 0.0 function opf_objective(x, param) #start = time() cost = 0.0 - for (i,gen) in ref[:gen] - pg = x[var_lookup["pg_$(i)"]] - cost += gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3] + for i in ref_gen_idxs + pg = x[lookup_pg[i]] + _cost_arr = cost_arrs[i] + cost += _cost_arr[1]*pg^2 + _cost_arr[2]*pg + _cost_arr[3] end #total_callback_time += time() - start return cost end function opf_constraints(ret, x, param) - #start = time() - va = Dict(i => x[var_lookup["va_$(i)"]] for (i,bus) in ref[:bus]) - vm = Dict(i => x[var_lookup["vm_$(i)"]] for (i,bus) in ref[:bus]) - - pg = Dict(i => x[var_lookup["pg_$(i)"]] for (i,gen) in ref[:gen]) - qg = Dict(i => x[var_lookup["qg_$(i)"]] for (i,gen) in ref[:gen]) - - p = Dict((l,i,j) => x[var_lookup["p_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs]) - q = Dict((l,i,j) => x[var_lookup["q_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs]) - - vm_fr = Dict(l => vm[branch["f_bus"]] for (l,branch) in ref[:branch]) - vm_to = Dict(l => vm[branch["t_bus"]] for (l,branch) in ref[:branch]) - va_fr = Dict(l => va[branch["f_bus"]] for (l,branch) in ref[:branch]) - va_to = Dict(l => va[branch["t_bus"]] for (l,branch) in ref[:branch]) - + offsetidx = 0 - va_con = [va[i] for (i,bus) in ref[:ref_buses]] + # va_con + for (reti,i) in enumerate(ref_buses_idxs) + ret[reti + offsetidx] = x[lookup_va[i]] + end + offsetidx += length(ref_buses_idxs) + # @constraint(model, # sum(p[a] for a in ref[:bus_arcs][i]) == - # sum(pg[g] for g in ref[:bus_gens][i]) - + # sum(pg[g] for g in ref_bus_gens[i]) - # sum(load["pd"] for load in bus_loads) - - # sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2 + # sum(shunt["gs"] for shunt in bus_shunts)*x[lookup_vm[i]]^2 # ) - power_balance_p_con = [ - sum(pg[j] for j in ref[:bus_gens][i]; init=0.0) - - bus_pd[i] - - bus_gs[i]*vm[i]^2 - - sum(p[a] for a in ref[:bus_arcs][i]) - for (i,bus) in ref[:bus] - ] + + # power_balance_p_con + for (reti,i) in enumerate(ref_bus_idxs) + ret[reti + offsetidx] = sum(x[lookup_pg[j]] for j in ref_bus_gens[i]; init=0.0) - + bus_pd[i] - + bus_gs[i]*x[lookup_vm[i]]^2 - + sum(x[p_idxmap[a]] for a in ref_bus_arcs[i]) + end + + offsetidx += length(ref_bus_idxs) # @constraint(model, # sum(q[a] for a in ref[:bus_arcs][i]) == - # sum(qg[g] for g in ref[:bus_gens][i]) - + # sum(x[lookup_qg[g]] for g in ref_bus_gens[i]) - # sum(load["qd"] for load in bus_loads) + - # sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2 + # sum(shunt["bs"] for shunt in bus_shunts)*x[lookup_vm[i]]^2 # ) - power_balance_q_con = [ - sum(qg[j] for j in ref[:bus_gens][i]; init=0.0) - - bus_qd[i] + - bus_bs[i]*vm[i]^2 - - sum(q[a] for a in ref[:bus_arcs][i]) - for (i,bus) in ref[:bus] - ] + # power_balance_q_con + for (reti,i) in enumerate(ref_bus_idxs) + ret[reti + offsetidx] = sum(x[lookup_qg[j]] for j in ref_bus_gens[i]; init=0.0) - + bus_qd[i] + + bus_bs[i]*x[lookup_vm[i]]^2 - + sum(x[q_idxmap[a]] for a in ref_bus_arcs[i]) + end + offsetidx += length(ref_bus_idxs) # @NLconstraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) ) - power_flow_p_from_con = [ - (br_g[l]+br_g_fr[l])/br_ttm[l]*vm_fr[l]^2 + - (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) + - (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) - - p[(l,i,j)] - for (l,i,j) in ref[:arcs_from] - ] + # power_flow_p_from_con = + for (reti,(l,i,j)) in enumerate(ref_arcs_from) + ret[reti + offsetidx] = (br_g[l]+br_g_fr[l])/br_ttm[l]*x[lookup_vm[f_bus[l]]]^2 + + (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*cos(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) + + (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*sin(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) - + x[p_idxmap[(l,i,j)]] + end + + offsetidx += length(ref_arcs_from) # @NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) ) - power_flow_p_to_con = [ - (br_g[l]+br_g_to[l])*vm_to[l]^2 + - (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) + - (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) - - p[(l,i,j)] - for (l,i,j) in ref[:arcs_to] - ] + # power_flow_p_to_con + for (reti,(l,i,j)) in enumerate(ref_arcs_to) + ret[reti + offsetidx] = (br_g[l]+br_g_to[l])*x[lookup_vm[t_bus[l]]]^2 + + (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*cos(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) + + (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*sin(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) - + x[p_idxmap[(l,i,j)]] + end + + offsetidx += length(ref_arcs_to) # @NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) ) - power_flow_q_from_con = [ - -(br_b[l]+br_b_fr[l])/br_ttm[l]*vm_fr[l]^2 - - (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) + - (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) - - q[(l,i,j)] - for (l,i,j) in ref[:arcs_from] - ] + # power_flow_q_from_con + for (reti,(l,i,j)) in enumerate(ref_arcs_from) + ret[reti + offsetidx] = -(br_b[l]+br_b_fr[l])/br_ttm[l]*x[lookup_vm[f_bus[l]]]^2 - + (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*cos(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) + + (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*sin(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) - + x[q_idxmap[(l,i,j)]] + end + + offsetidx += length(ref_arcs_from) # @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) ) - power_flow_q_to_con = [ - -(br_b[l]+br_b_to[l])*vm_to[l]^2 - - (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) + - (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) - - q[(l,i,j)] - for (l,i,j) in ref[:arcs_to] - ] + # power_flow_q_to_con + for (reti,(l,i,j)) in enumerate(ref_arcs_to) + ret[reti + offsetidx] = -(br_b[l]+br_b_to[l])*x[lookup_vm[t_bus[l]]]^2 - + (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*cos(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) + + (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*sin(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) - + x[q_idxmap[(l,i,j)]] + end + + offsetidx += length(ref_arcs_to) # @constraint(model, va_fr - va_to <= branch["angmax"]) # @constraint(model, va_fr - va_to >= branch["angmin"]) - power_flow_vad_con = [ - va_fr[l] - va_to[l] - for (l,i,j) in ref[:arcs_from] - ] + # power_flow_vad_con + for (reti,(l,i,j)) in enumerate(ref_arcs_from) + ret[reti + offsetidx] = x[lookup_va[f_bus[l]]] - x[lookup_va[t_bus[l]]] + end + + offsetidx += length(ref_arcs_from) # @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2) - power_flow_mva_from_con = [ - p[(l,i,j)]^2 + q[(l,i,j)]^2 - for (l,i,j) in ref[:arcs_from] - ] + # power_flow_mva_from_con + for (reti,(l,i,j)) in enumerate(ref_arcs_from) + ret[reti + offsetidx] = x[p_idxmap[(l,i,j)]]^2 + x[q_idxmap[(l,i,j)]]^2 + end + + offsetidx += length(ref_arcs_from) # @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2) - power_flow_mva_to_con = [ - p[(l,i,j)]^2 + q[(l,i,j)]^2 - for (l,i,j) in ref[:arcs_to] - ] - ret .= [ - va_con..., - power_balance_p_con..., - power_balance_q_con..., - power_flow_p_from_con..., - power_flow_p_to_con..., - power_flow_q_from_con..., - power_flow_q_to_con..., - power_flow_vad_con..., - power_flow_mva_from_con..., - power_flow_mva_to_con..., - ] - #total_callback_time += time() - start + # power_flow_mva_to_con + for (reti,(l,i,j)) in enumerate(ref_arcs_to) + ret[reti + offsetidx] = x[p_idxmap[(l,i,j)]]^2 + x[q_idxmap[(l,i,j)]]^2 + end + + offsetidx += length(ref_arcs_to) + + @assert offsetidx == length(ret) + nothing end con_lbs = Float64[] @@ -331,56 +521,60 @@ function solve_opf_optimization(file_name) model_variables = length(var_init) ret = Array{Float64}(undef, length(con_lbs)) - model_constraints = length(opf_constraints(ret, var_init, ref)) - println("variables: $(model_variables), $(length(var_lb)), $(length(var_ub))") - println("constraints: $(model_constraints), $(length(con_lbs)), $(length(con_ubs))") - + model_constraints = length(con_lbs) - optprob = Optimization.OptimizationFunction(opf_objective, Optimization.AutoModelingToolkit(true, true); cons=opf_constraints) + optprob = Optimization.OptimizationFunction(opf_objective, adchoice; cons=opf_constraints) prob = Optimization.OptimizationProblem(optprob, var_init; lb=var_lb, ub=var_ub, lcons=con_lbs, ucons=con_ubs) +end - model_build_time = time() - time_model_start - +function solve_opf_optimization(dataset; adchoice = Optimization.AutoForwardDiff()) + model_build_time = @elapsed prob = build_opf_optimization_prob(dataset; adchoice) - time_solve_start = time() + # Correctness tests + @assert prob.f(prob.u0, nothing) == 0.0 + ret = zeros(length(prob.lcons)) + prob.f.cons(ret, prob.u0, nothing) + @allocated prob.f(prob.u0, nothing) == 0 + @allocated prob.f.cons(ret, prob.u0, nothing) == 0 - sol = Optimization.solve(prob, Ipopt.Optimizer()) + solve_time_with_compilation = @elapsed sol = Optimization.solve(prob, Ipopt.Optimizer()) cost = sol.minimum feasible = (sol.retcode == Optimization.SciMLBase.ReturnCode.Success) #println(sol.u) # solution vector - solve_time = time() - time_solve_start - total_time = time() - time_data_start - - println("") - println("\033[1mSummary\033[0m") - println(" case........: $(file_name)") - println(" variables...: $(model_variables)") - println(" constraints.: $(model_constraints)") - println(" feasible....: $(feasible)") - println(" cost........: $(round(Int, cost))") - println(" total time..: $(total_time)") - println(" data time.: $(data_load_time)") - println(" build time: $(model_build_time)") - println(" solve time: $(solve_time)") - # println(" callbacks: $(total_callback_time)") - println("") - - return Dict( + solve_time_without_compilation = @elapsed sol = Optimization.solve(prob, Ipopt.Optimizer()) + + return (prob,sol),Dict( "case" => file_name, - "variables" => model_variables, - "constraints" => model_constraints, + "variables" => length(prob.u0), + "constraints" => length(prob.lcons), "feasible" => feasible, "cost" => cost, - "time_total" => total_time, - "time_data" => data_load_time, "time_build" => model_build_time, - "time_solve" => solve_time, - #"time_callbacks" => TBD, + "time_solve" => solve_time_without_compilation, + "time_solve_compilation" => solve_time_with_compilation, ) end -solve_opf_optimization("../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m") +function test_optimization_prob(dataset, test_u0) + prob = build_opf_optimization_prob(dataset) + ret = zeros(length(prob.lcons)) + prob.f.cons(ret, test_u0, nothing) + obj = prob.f(test_u0, nothing) + obj, ret +end +``` + +```julia +optimization_test_res = test_optimization_prob(dataset, test_u0) +``` + +```julia +@assert optimization_test_res[1] == test_obj +``` + +```julia +@assert optimization_test_res[2] == test_cons ``` ## JuMP.jl @@ -393,47 +587,42 @@ import PowerModels import Ipopt import JuMP -function solve_opf_jump(file_name) - - - - time_model_start = time() - +function build_opf_jump_prob(dataset) + (;data, ref) = dataset + constraints = Any[] model = JuMP.Model(Ipopt.Optimizer) #JuMP.set_optimizer_attribute(model, "print_level", 0) - JuMP.@variable(model, va[i in keys(ref[:bus])]) - JuMP.@variable(model, ref[:bus][i]["vmin"] <= vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0) - - JuMP.@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"]) - JuMP.@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"]) - - JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"]) - JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"]) + vars = [JuMP.@variable(model, va[i in keys(ref[:bus])]), + JuMP.@variable(model, ref[:bus][i]["vmin"] <= vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0), + JuMP.@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"]), + JuMP.@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"]), + JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"]), + JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])] JuMP.@objective(model, Min, sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen])) for (i,bus) in ref[:ref_buses] - JuMP.@constraint(model, va[i] == 0) + push!(constraints,JuMP.@constraint(model, va[i] == 0)) end for (i,bus) in ref[:bus] bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]] bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]] - JuMP.@constraint(model, + push!(constraints,JuMP.@constraint(model, sum(p[a] for a in ref[:bus_arcs][i]) == sum(pg[g] for g in ref[:bus_gens][i]) - sum(load["pd"] for load in bus_loads) - sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2 - ) + )) - JuMP.@constraint(model, + push!(constraints,JuMP.@constraint(model, sum(q[a] for a in ref[:bus_arcs][i]) == sum(qg[g] for g in ref[:bus_gens][i]) - sum(load["qd"] for load in bus_loads) + sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2 - ) + )) end # Branch power flow physics and limit constraints @@ -460,19 +649,19 @@ function solve_opf_jump(file_name) b_to = branch["b_to"] # From side of the branch flow - JuMP.@NLconstraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) ) - JuMP.@NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) ) + push!(constraints,JuMP.@NLconstraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )) + push!(constraints,JuMP.@NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )) # To side of the branch flow - JuMP.@NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) ) - JuMP.@NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) ) + push!(constraints,JuMP.@NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )) + push!(constraints,JuMP.@NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )) # Voltage angle difference limit - JuMP.@constraint(model, branch["angmin"] <= va_fr - va_to <= branch["angmax"]) + push!(constraints,JuMP.@constraint(model, branch["angmin"] <= va_fr - va_to <= branch["angmax"])) # Apparent power limit, from side and to side - JuMP.@constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2) - JuMP.@constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2) + push!(constraints,JuMP.@constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)) + push!(constraints,JuMP.@constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)) end model_variables = JuMP.num_variables(model) @@ -481,18 +670,17 @@ function solve_opf_jump(file_name) non_nl_constraints = sum(JuMP.num_constraints(model, ft, st) for (ft, st) in JuMP.list_of_constraint_types(model) if ft != JuMP.VariableRef) model_constraints = JuMP.num_nonlinear_constraints(model) + non_nl_constraints - model_build_time = time() - time_model_start - + model, vars, constraints +end - time_solve_start = time() +function solve_opf_jump(dataset) + model_build_time = @elapsed model = build_opf_jump_prob(dataset)[1] + solve_time_with_compilation = @elapsed JuMP.optimize!(model) + solve_time_without_compilation = @elapsed JuMP.optimize!(model) - JuMP.optimize!(model) cost = JuMP.objective_value(model) feasible = (JuMP.termination_status(model) == JuMP.LOCALLY_SOLVED) - solve_time = time() - time_solve_start - total_time = time() - time_data_start - nlp_block = JuMP.MOI.get(model, JuMP.MOI.NLPBlock()) total_callback_time = nlp_block.evaluator.eval_objective_timer + @@ -500,43 +688,87 @@ function solve_opf_jump(file_name) nlp_block.evaluator.eval_constraint_timer + nlp_block.evaluator.eval_constraint_jacobian_timer + nlp_block.evaluator.eval_hessian_lagrangian_timer - - println("") - println("\033[1mSummary\033[0m") - println(" case........: $(file_name)") - println(" variables...: $(model_variables)") - println(" constraints.: $(model_constraints)") - println(" feasible....: $(feasible)") - println(" cost........: $(round(Int, cost))") - println(" total time..: $(total_time)") - println(" data time.: $(data_load_time)") - println(" build time: $(model_build_time)") - println(" solve time: $(solve_time)") - println(" callbacks: $(total_callback_time)") - println("") - println(" callbacks time:") - println(" * obj.....: $(nlp_block.evaluator.eval_objective_timer)") - println(" * grad....: $(nlp_block.evaluator.eval_objective_gradient_timer)") - println(" * cons....: $(nlp_block.evaluator.eval_constraint_timer)") - println(" * jac.....: $(nlp_block.evaluator.eval_constraint_jacobian_timer)") - println(" * hesslag.: $(nlp_block.evaluator.eval_hessian_lagrangian_timer)") - println("") - - return Dict( + model_variables = JuMP.num_variables(model) + non_nl_constraints = sum(JuMP.num_constraints(model, ft, st) for (ft, st) in JuMP.list_of_constraint_types(model) if ft != JuMP.VariableRef) + model_constraints = JuMP.num_nonlinear_constraints(model) + non_nl_constraints + + return model, Dict( "case" => file_name, "variables" => model_variables, "constraints" => model_constraints, "feasible" => feasible, "cost" => cost, - "time_total" => total_time, - "time_data" => data_load_time, "time_build" => model_build_time, - "time_solve" => solve_time, - "time_callbacks" => total_callback_time, + "time_solve" => solve_time_without_compilation, + "time_solve_compilation" => solve_time_with_compilation, ) end -solve_opf_jump("../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m") +function test_jump_prob(dataset, test_u0) + model, vars, constraints = build_opf_jump_prob(dataset) + (; + lookup_pg, + lookup_qg, + lookup_va, + lookup_vm, + lookup_lij, + lookup_p_lij, + lookup_q_lij) = dataset + f = JuMP.objective_function(model) + + flatvars = reduce(vcat,[reduce(vcat,vars[i]) for i in 1:length(vars)]) + point = Dict() + for v in flatvars + varname, varint = split(JuMP.name(v), "[") + idx = if varint[1] == '(' + varint = (parse(Int, varint[2]), parse(Int, varint[5]), parse(Int, varint[8])) + if varname == "p" + lookup_p_lij[findfirst(x->x==varint,lookup_lij)] + elseif varname == "q" + lookup_q_lij[findfirst(x->x==varint,lookup_lij)] + else + error("Invalid $varname, $varint") + end + else + varint = parse(Int, varint[1:end-1]) + if varname == "va" + lookup_va[varint] + elseif varname == "pg" + lookup_pg[varint] + elseif varname == "qg" + lookup_qg[varint] + elseif varname == "vm" + lookup_vm[varint] + else + error("Invalid $varname, $varint") + end + end + point[v] = test_u0[idx] + end + obj = JuMP.value(x->point[x], f) + cons = [JuMP.value(x->point[x], c) for c in constraints] + obj, cons +end +``` + +```julia +jump_test_res = test_jump_prob(dataset, test_u0) +``` + +```julia +@assert jump_test_res[1] ≈ test_obj +``` + +```julia +@assert sort(abs.(jump_test_res[2])) ≈ sort(abs.(test_cons)) +``` + +```julia +println(sort(abs.(jump_test_res[2]))) +``` + +```julia +println(sort(abs.(test_cons))) ``` ## NLPModels.jl @@ -545,17 +777,11 @@ Implementation reference: https://juliasmoothoptimizers.github.io/ADNLPModels.jl Other AD libraries can be considered: https://juliasmoothoptimizers.github.io/ADNLPModels.jl/stable/ ```julia -import PowerModels -import Symbolics import ADNLPModels import NLPModelsIpopt -function solve_opf_nlpmodels(file_name) - - - - time_model_start = time() - +function build_opf_nlpmodels_prob(dataset) + (;data, ref) = dataset bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus]) bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus]) @@ -844,53 +1070,50 @@ function solve_opf_nlpmodels(file_name) model_variables = length(var_init) model_constraints = length(opf_constraints!(similar(con_lbs), var_init)) - println("variables: $(model_variables), $(length(var_lb)), $(length(var_ub))") - println("constraints: $(model_constraints), $(length(con_lbs)), $(length(con_ubs))") nlp = ADNLPModels.ADNLPModel!(opf_objective, var_init, var_lb, var_ub, opf_constraints!, con_lbs, con_ubs, backend = :optimized) +end - model_build_time = time() - time_model_start - - - time_solve_start = time() - - output = NLPModelsIpopt.ipopt(nlp) +function solve_opf_nlpmodels(dataset) + model_build_time = @elapsed nlp = build_opf_nlpmodels_prob(dataset) + solve_time = @elapsed output = NLPModelsIpopt.ipopt(nlp) cost = output.objective - feasible = (output.primal_feas <= 1e-6) - - solve_time = time() - time_solve_start - total_time = time() - time_data_start - + feasible = (output.primal_feas <= 1e-6) - println("") - println("\033[1mSummary\033[0m") - println(" case........: $(file_name)") - println(" variables...: $(model_variables)") - println(" constraints.: $(model_constraints)") - println(" feasible....: $(feasible)") - println(" cost........: $(round(Int, cost))") - println(" total time..: $(total_time)") - println(" data time.: $(data_load_time)") - println(" build time: $(model_build_time)") - println(" solve time: $(solve_time)") - # println(" callbacks: $(total_callback_time)") - println("") + model_variables = nlp.meta.nvar + model_constraints = nlp.meta.ncon - return Dict( + return (nlp, output), Dict( "case" => file_name, "variables" => model_variables, "constraints" => model_constraints, "feasible" => feasible, "cost" => cost, - "time_total" => total_time, - "time_data" => data_load_time, "time_build" => model_build_time, "time_solve" => solve_time, #"time_callbacks" => TBD, ) end -solve_opf_nlpmodels("../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m") +function test_nlpmodels_prob(dataset, test_u0) + nlp = build_opf_nlpmodels_prob(dataset) + ret = zeros(nlp.meta.ncon) + nlp.c!(ret, test_u0) + obj = nlp.f(test_u0) + obj, ret +end +``` + +```julia +nlpmodels_test_res = test_nlpmodels_prob(dataset, test_u0) +``` + +```julia +@assert nlpmodels_test_res[1] == test_obj +``` + +```julia +@assert nlpmodels_test_res[2] == test_cons ``` ## Nonconvex @@ -899,15 +1122,11 @@ Implementation reference: https://julianonconvex.github.io/Nonconvex.jl/stable/p Currently does not converge due to an upstream issue with the AD backend Zygote: https://github.com/JuliaNonconvex/Nonconvex.jl/issues/130 ```julia -import PowerModels import Nonconvex Nonconvex.@load Ipopt - -function solve_opf_nonconvex(file_name) - - - +function build_opf_nonconvex_prob(dataset) + (;data, ref) = dataset time_model_start = time() model = Nonconvex.DictModel() @@ -1100,18 +1319,13 @@ function solve_opf_nonconvex(file_name) add_ineq_constraint!(model, x -> const_voltage_angle_difference_lb(x,l,i,j)) add_ineq_constraint!(model, x -> const_voltage_angle_difference_ub(x,l,i,j)) end + model +end - model_variables = Nonconvex.NonconvexCore.getnvars(model) - model_constraints = Nonconvex.NonconvexCore.getnconstraints(model) - println("variables: $(model_variables)") - println("constraints: $(model_constraints)") - - model_build_time = time() - time_model_start - - - time_solve_start = time() +function solve_opf_nonconvex(dataset) + model_build_time = @elapsed model = build_opf_nonconvex_prob(dataset) - result = Nonconvex.optimize( + solve_time = @elapsed result = Nonconvex.optimize( model, IpoptAlg(), NonconvexCore.getinit(model); @@ -1121,38 +1335,80 @@ function solve_opf_nonconvex(file_name) cost = result.minimum feasible = result.status == 0 # just guessing this is correct for Ipopt - solve_time = time() - time_solve_start - total_time = time() - time_data_start - - println("") - println("\033[1mSummary\033[0m") - println(" case........: $(file_name)") - println(" variables...: $(model_variables)") - println(" constraints.: $(model_constraints)") - println(" feasible....: $(feasible)") - println(" cost........: $(round(Int, cost))") - println(" total time..: $(total_time)") - println(" data time.: $(data_load_time)") - println(" build time: $(model_build_time)") - println(" solve time: $(solve_time)") - # println(" callbacks: $(total_callback_time)") - println("") + model_variables = Nonconvex.NonconvexCore.getnvars(model) + model_constraints = Nonconvex.NonconvexCore.getnconstraints(model) - return Dict( + return (model, result), Dict( "case" => file_name, "variables" => model_variables, "constraints" => model_constraints, "feasible" => feasible, "cost" => cost, - "time_total" => total_time, - "time_data" => data_load_time, "time_build" => model_build_time, "time_solve" => solve_time, #"time_callbacks" => TBD, ) end -solve_opf_nonconvex("../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m") +function test_nonconvex_prob(dataset, test_u0) + model = build_opf_nonconvex_prob(dataset) + (; + lookup_pg, + lookup_qg, + lookup_va, + lookup_vm, + lookup_lij, + lookup_p_lij, + lookup_q_lij) = dataset + + point = Dict() + for v in keys(model.init) + varsplit = split(v, "_") + varname = varsplit[1] + varint = parse.(Int, varsplit[2:end]) + + idx = if varname == "p" + lookup_p_lij[findfirst(x->x==Tuple(varint),lookup_lij)] + elseif varname == "q" + lookup_q_lij[findfirst(x->x==Tuple(varint),lookup_lij)] + elseif varname == "va" + lookup_va[varint[1]] + elseif varname == "pg" + lookup_pg[varint[1]] + elseif varname == "qg" + lookup_qg[varint[1]] + elseif varname == "vm" + lookup_vm[varint[1]] + else + error("Invalid $varname, $varint") + end + point[v] = test_u0[idx] + end + u0 = OrderedDict(keys(model.init) .=> getindex.((point,),keys(model.init))) + obj = model.objective(u0) + cons = vcat(model.eq_constraints(u0), model.ineq_constraints(u0)) + obj, cons +end +``` + +```julia +nonconvex_test_res = test_nonconvex_prob(dataset, test_u0) +``` + +```julia +@assert nonconvex_test_res[1] ≈ test_obj +``` + +```julia +@assert sort(abs.(nonconvex_test_res[2])) ≈ sort(abs.(test_cons)) +``` + +```julia +println(sort(abs.(nonconvex_test_res[2]))) +``` + +```julia +println(sort(abs.(test_cons))) ``` ## Optim.jl @@ -1162,14 +1418,10 @@ Currently does not converge to a feasible point, root cause in unclear `debug/optim-debug.jl` can be used to confirm it will converge if given a suitable starting point ```julia -import PowerModels import Optim -function solve_opf_optim(file_name) - - - - time_model_start = time() +function build_opf_optim_prob(dataset) + (;data, ref) = dataset bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus]) bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus]) @@ -1476,78 +1728,58 @@ function solve_opf_optim(file_name) push!(con_ubs, branch["rate_a"]^2) end - model_variables = length(var_init) - model_constraints = length(opf_constraints(zeros(length(con_lbs)), var_init)) - println("variables: $(model_variables), $(length(var_lb)), $(length(var_ub))") - println("constraints: $(model_constraints), $(length(con_lbs)), $(length(con_ubs))") - df = Optim.TwiceDifferentiable(opf_objective, var_init) dfc = Optim.TwiceDifferentiableConstraints(opf_constraints, var_lb, var_ub, con_lbs, con_ubs) + df, dfc, var_init, con_lbs, con_ubs +end - model_build_time = time() - time_model_start - - time_solve_start = time() +function solve_opf_optim(dataset) + model_build_time = @elapsed df, dfc, var_init, con_lbs, con_ubs = build_opf_optim_prob(dataset) options = Optim.Options(show_trace=true) - - # NOTE: had to change initial guess to be an interior point, otherwise getting NaN values - res = Optim.optimize(df, dfc, var_init, Optim.IPNewton(), options) - #res = Optim.optimize(df, dfc, var_init, Optim.LBFGS(), options) # StackOverflowError: - #res = Optim.optimize(df, dfc, var_init, Optim.NelderMead(), options) # StackOverflowError: - display(res) + solve_time = @elapsed res = Optim.optimize(df, dfc, var_init, Optim.IPNewton(), options) sol = res.minimizer cost = res.minimum - solve_time = time() - time_solve_start - total_time = time() - time_data_start - - # NOTE: confirmed these constraint violations can be eliminated # if a better starting point is used - sol_eval = opf_constraints(zeros(length(con_lbs)), sol) + sol_eval = dfc.c!(zeros(dfc.bounds.nc), sol) vio_lb = [max(v,0) for v in (con_lbs .- sol_eval)] vio_ub = [max(v,0) for v in (sol_eval .- con_ubs)] const_vio = vio_lb .+ vio_ub - #println(const_vio) - println("total constraint violation: $(sum(const_vio))") constraint_tol = 1e-6 feasible = (sum(const_vio) <= constraint_tol) - if !feasible - @warn "Optim optimize failed to satify the problem constraints" - end - - println("") - println("\033[1mSummary\033[0m") - println(" case........: $(file_name)") - println(" variables...: $(model_variables)") - println(" constraints.: $(model_constraints)") - println(" feasible....: $(feasible)") - println(" cost........: $(round(Int, cost))") - println(" total time..: $(total_time)") - println(" data time.: $(data_load_time)") - println(" build time: $(model_build_time)") - println(" solve time: $(solve_time)") - # println(" callbacks: $(total_callback_time)") - println("") - - - return Dict( + return (res,), Dict( "case" => file_name, - "variables" => model_variables, - "constraints" => model_constraints, + "variables" => length(var_init), + "constraints" => dfc.bounds.nc, "feasible" => feasible, "cost" => cost, - "time_total" => total_time, - "time_data" => data_load_time, "time_build" => model_build_time, "time_solve" => solve_time, - #"time_callbacks" => TBD, ) end -solve_opf_optim("../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m") +function test_optim_prob(dataset, test_u0) + df, dfc, var_init, con_lbs, con_ubs = build_opf_optim_prob(dataset) + obj = df.f(test_u0) + cons = dfc.c!(zeros(dfc.bounds.nc), test_u0) + obj, cons +end +``` + +```julia +optim_test_res = test_optim_prob(dataset, test_u0) +``` + +```julia +@assert optim_test_res[1] == test_obj +``` + +```julia +@assert optim_test_res[2] == test_cons ``` ## CASADI @@ -1562,8 +1794,8 @@ import PythonCall import CondaPkg CondaPkg.add("casadi") -function solve_opf_casadi(file_name) - +function solve_opf_casadi(dataset) + (;data, ref) = dataset time_model_start = time() @@ -1737,4 +1969,140 @@ function solve_opf_casadi(file_name) end solve_opf_casadi("../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m") +``` + +## Start the Benchmarking + +```julia +file_name = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m" +dataset = load_and_setup_data(file_name); +``` + +```julia +model, res = solve_opf_optimization(dataset); +res +``` + +```julia +model, res = solve_opf_jump(dataset); +res +``` + +```julia +model, res = solve_opf_nlpmodels(dataset); +res +``` + +```julia +model, res = solve_opf_nonconvex(dataset); +res +``` + +```julia +model, res = solve_opf_optim(dataset); +res +``` + +```julia +file_name = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m" +dataset = load_and_setup_data(file_name); +``` + +```julia +model, res = solve_opf_optimization(dataset); +res +``` + +```julia +model, res = solve_opf_jump(dataset); +res +``` + +```julia +model, res = solve_opf_nlpmodels(dataset); +res +``` + +```julia +model, res = solve_opf_nonconvex(dataset); +res +``` + +```julia +model, res = solve_opf_optim(dataset); +res +``` + +```julia +using DataFrames, PrettyTables + +function multidata_multisolver_benchmark(dataset_files) + + cases = String[] + vars = Int[] + cons = Int[] + optimization_time = Float64[] + jump_time = Float64[] + nlpmodels_time = Float64[] + nonconvex_time = Float64[] + optim_time = Float64[] + + optimization_cost = Float64[] + jump_cost = Float64[] + nlpmodels_cost = Float64[] + nonconvex_cost = Float64[] + optim_cost = Float64[] + + for file in dataset_files + @show file + dataset = load_and_setup_data(file) + model, res = solve_opf_optimization(dataset) + push!(cases, split(res["case"],"/")[end]) + push!(vars, res["variables"]) + push!(cons, res["constraints"]) + push!(optimization_time, res["time_solve"]) + push!(optimization_cost, res["cost"]) + + model, res = solve_opf_jump(dataset) + push!(jump_time, res["time_solve"]) + push!(jump_cost, res["cost"]) + + model, res = solve_opf_nlpmodels(dataset) + push!(nlpmodels_time, res["time_solve"]) + push!(nlpmodels_cost, res["cost"]) + + model, res = solve_opf_nonconvex(dataset) + push!(nonconvex_time, res["time_solve"]) + push!(nonconvex_cost, res["cost"]) + + model, res = solve_opf_optim(dataset) + push!(optim_time, res["time_solve"]) + push!(optim_cost, res["cost"]) + end + DataFrame(:case => cases, :vars => vars, :cons => cons, + :optimization => optimization_time, :optimization_cost => optimization_cost, + :jump => jump_time, :jump_cost => jump_cost, + :nlpmodels => nlpmodels_time, :nlpmodels_cost => nlpmodels_cost, + :nonconvex => nonconvex_time, :nonconvex_cost => nonconvex_cost, + :optim => optim_time, :optim_cost => optim_cost) +end + +test_datasets = [ + "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m", + "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m" + ] +``` + +```julia +timing_data = multidata_multisolver_benchmark(test_datasets) +``` + +```julia +pretty_table(timing_data) +``` + +```julia +#file_name = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case118_ieee.m" +#dataset = load_and_setup_data(file_name); +#build_and_solve_optimization(dataset) ``` \ No newline at end of file