[Do not merge] a script to benchmark impact of lazy PTDF

Project.toml

@@ -13,8 +13,10 @@ JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+MathOptLazy = "5d5fe9b5-b0a4-4485-81f6-7b1b939155e1"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -23,14 +25,15 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 [weakdeps]
 PowerFlows = "94fada2c-0ca5-4b90-a1fb-4bc5b59ccfc7"
 
+[sources]
+InfrastructureOptimizationModels = {rev = "lk/pom-test-fixes", url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl"}
+InfrastructureSystems = {rev = "IS4", url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl"}
+MathOptLazy = {rev = "main", url = "https://github.com/jump-dev/MathOptLazy.jl"}
+PowerSystems = {rev = "psy6", url = "https://github.com/NREL-Sienna/PowerSystems.jl"}
+
 [extensions]
 PowerFlowsExt = "PowerFlows"
 
-[sources]
-InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"}
-PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "psy6"}
-InfrastructureOptimizationModels = {url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl", rev = "lk/pom-test-fixes"}
-
 [compat]
 Dates = "1"
 DocStringExtensions = "~0.8, ~0.9"
@@ -40,6 +43,7 @@ InteractiveUtils = "1.11.0"
 JuMP = "^1.28"
 PowerNetworkMatrices = "^0.19"
 PowerSystems = "5.3"
+PrettyTables = "3"
 ProgressMeter = "1.11.0"
 TimerOutputs = "~0.5"
 julia = "^1.11"

src/PowerModels/core/objective.jl

@@ -157,7 +157,7 @@ function expression_pg_cost(pm::AbstractPowerModel; report::Bool = true)
                 )
             else
                 error(
-                    "Only cost models of types 1 and 2 are supported at this time, given cost model type of $(model) on generator $(i)",
+                    "Only cost models of types 1 and 2 are supported at this time, given cost model type of $(gen["model"]) on generator $(i)",
                 )
             end
         end

src/PowerOperationsModels.jl

@@ -9,8 +9,10 @@ import InfrastructureSystems: @assert_op, TableFormat
 import JuMP
 import JuMP.Containers: DenseAxisArray, SparseAxisArray
 import Logging
+import MathOptLazy
 import PowerNetworkMatrices
 import ProgressMeter
+import PrettyTables
 import PowerSystems
 import PowerSystems: get_component
 import Serialization

src/ac_transmission_models/AC_branches.jl

@@ -362,6 +362,21 @@ function get_min_max_limits(
     return (min = -π / 2, max = π / 2)
 end
 
+_get_tag(model::JuMP.GenericModel) = _get_tag(JuMP.backend(model))
+
+_get_tag(::JuMP.MOI.ModelLike) = ()
+
+_get_tag(::MathOptLazy.Optimizer) = (MathOptLazy.Lazy(),)
+
+_get_tag(model::JuMP.MOI.Bridges.LazyBridgeOptimizer) = _get_tag(model.model)
+
+function _get_tag(model::JuMP.MOI.Utilities.CachingOptimizer)
+    if JuMP.MOI.Utilities.state(model) == JuMP.MOI.Utilities.NO_OPTIMIZER
+        return ()
+    end
+    return _get_tag(model.optimizer)
+end
+
 function _add_flow_rate_constraint!(
     container::OptimizationContainer,
     ::Type{T},
@@ -373,24 +388,32 @@ function _add_flow_rate_constraint!(
     branch_maps_by_type::Dict,
     name::String,
 ) where {T <: PSY.ACTransmission}
-    reduction_entry = branch_maps_by_type[arc]
     time_steps = get_time_steps(container)
+    limits = get_min_max_limits(branch_maps_by_type[arc], FlowRateConstraint, StaticBranch)
+    model = get_jump_model(container)
+    tag = _get_tag(model)
     if use_slacks
-        slack_ub = get_variable(container, FlowActivePowerSlackUpperBound, T)[name, :]
-        slack_lb = get_variable(container, FlowActivePowerSlackLowerBound, T)[name, :]
-    end
-    limits = get_min_max_limits(reduction_entry, FlowRateConstraint, StaticBranch)
-    for t in time_steps
-        con_ub[name, t] =
-            JuMP.@constraint(
-                get_jump_model(container),
-                var[name, t] - (use_slacks ? slack_ub[t] : 0.0) <= limits.max
-            )
-        con_lb[name, t] =
-            JuMP.@constraint(
-                get_jump_model(container),
-                var[name, t] + (use_slacks ? slack_lb[t] : 0.0) >= limits.min
-            )
+        slack_ub = get_variable(container, FlowActivePowerSlackUpperBound, T)
+        slack_lb = get_variable(container, FlowActivePowerSlackLowerBound, T)
+        for t in time_steps
+            con_ub[name, t] =
+                JuMP.@constraint(
+                    model,
+                    var[name, t] - slack_ub[name, t] <= limits.max,
+                    tag...
+                )
+            con_lb[name, t] =
+                JuMP.@constraint(
+                    model,
+                    var[name, t] + slack_lb[name, t] >= limits.min,
+                    tag...
+                )
+        end
+    else
+        for t in time_steps
+            con_ub[name, t] = JuMP.@constraint(model, var[name, t] <= limits.max, tag...)
+            con_lb[name, t] = JuMP.@constraint(model, var[name, t] >= limits.min, tag...)
+        end
     end
     return
 end

src/network_models/instantiate_network_model.jl

@@ -75,7 +75,7 @@ function _get_irreducible_buses_due_to_dlrs(
     irreducible_buses = Set{Int64}()
     for branch_type in network_model.modeled_branch_types
         branch_type <: PSY.ACTransmission || continue
-        device_model = branch_models[Symbol(branch_type)]
+        device_model = branch_models[nameof(branch_type)]
         if !haskey(
             get_time_series_names(device_model),
             DynamicBranchRatingTimeSeriesParameter,

src/operation/decision_model.jl

@@ -67,7 +67,7 @@ function build!(
     IOM.register_recorders!(model, file_mode)
     logger = IS.configure_logging(get_internal(model), IOM.PROBLEM_LOG_FILENAME, file_mode)
     if store_system_in_results
-        @warn "store_system_in_results for $(model) is set to true. This will do nothing unless a Simulation is being built."
+        @warn "store_system_in_results is set to true. This will do nothing unless a Simulation is being built."
     end
     try
         Logging.with_logger(logger) do
@@ -151,7 +151,7 @@ function solve!(
     kwargs...,
 )
     if store_system_in_results
-        @warn "store_system_in_results for $(model) is set to true. This will do nothing unless a Simulation is being built."
+        @warn "store_system_in_results is set to true. This will do nothing unless a Simulation is being built."
     end
     build_if_not_already_built!(
         model;

src/operation/emulation_model.jl

@@ -69,7 +69,7 @@ function build!(
         file_mode,
     )
     if store_system_in_results
-        @warn "store_system_in_results for $(model) is set to true. This will do nothing unless a Simulation is being built."
+        @warn "store_system_in_results is set to true. This will do nothing unless a Simulation is being built."
     end
     try
         Logging.with_logger(logger) do
@@ -197,7 +197,7 @@ function run!(
     kwargs...,
 )
     if store_system_in_results
-        @warn "store_system_in_results for $(model) is set to true. This will do nothing unless a Simulation is being built."
+        @warn "store_system_in_results is set to true. This will do nothing unless a Simulation is being built."
     end
     build_if_not_already_built!(
         model;

test/Project.toml

@@ -5,16 +5,19 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 InfrastructureOptimizationModels = "bed98974-b02a-5e2f-9ee0-a103f5c45069"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+MathOptLazy = "5d5fe9b5-b0a4-4485-81f6-7b1b939155e1"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
@@ -32,10 +35,11 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [sources]
-InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"}
-PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "psy6"}
-InfrastructureOptimizationModels = {url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl", rev = "lk/pom-test-fixes"}
-PowerSystemCaseBuilder = {url = "https://github.com/NREL-Sienna/PowerSystemCaseBuilder.jl", rev = "psy6"}
+InfrastructureOptimizationModels = {rev = "lk/pom-test-fixes", url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl"}
+InfrastructureSystems = {rev = "IS4", url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl"}
+MathOptLazy = {rev = "main", url = "https://github.com/jump-dev/MathOptLazy.jl"}
+PowerSystemCaseBuilder = {rev = "psy6", url = "https://github.com/NREL-Sienna/PowerSystemCaseBuilder.jl"}
+PowerSystems = {rev = "psy6", url = "https://github.com/NREL-Sienna/PowerSystems.jl"}
 
 [compat]
 HiGHS = "1"

test/performance/simple_cats.jl

@@ -0,0 +1,147 @@
+using Revise
+import Dates
+import Gurobi
+import MathOptLazy
+import PowerNetworkMatrices as PNM
+import PowerOperationsModels as POM
+import PowerSystems as PSY
+import HiGHS
+using JuMP
+using PowerSystems
+
+function run_problem(optimizer)
+    sys = PSY.System(joinpath(@__DIR__, "CATS_Sienna.json"))
+    PSY.transform_single_time_series!(sys, Dates.Hour(24), Dates.Hour(24))
+    ptdf = POM.VirtualPTDF(
+        sys;
+        tol = 0.01,
+        network_reductions = [PNM.RadialReduction(), PNM.DegreeTwoReduction()],
+    )
+    network_model = POM.NetworkModel(
+        POM.PTDFPowerModel;
+        PTDF_matrix = ptdf,
+        reduce_radial_branches = true,
+        reduce_degree_two_branches = true,
+    )
+    template = POM.OperationsProblemTemplate(network_model)
+    POM.set_device_model!(template, PSY.ThermalStandard, POM.ThermalBasicUnitCommitment)
+    POM.set_device_model!(template, PSY.RenewableDispatch, POM.RenewableFullDispatch)
+    POM.set_device_model!(template, PSY.HydroDispatch, POM.HydroDispatchRunOfRiver)
+    POM.set_device_model!(template, PSY.PowerLoad, POM.StaticPowerLoad)
+    POM.set_device_model!(
+        template,
+        POM.DeviceModel(PSY.Line, POM.StaticBranch; use_slacks = true),
+    )
+    POM.set_device_model!(template, PSY.Transformer2W, POM.StaticBranch)
+    model = POM.DecisionModel(
+        template,
+        sys;
+        name = "CATS_UC2",
+        optimizer,
+        direct_mode_optimizer = true,
+        optimizer_solve_log_print = true,
+    )
+    POM.build!(model; output_dir = mktempdir(; cleanup = true))
+    @time POM.solve!(model)
+    return
+end
+
+# | Solver | Lazy? | tol=1e-2 | tol=1e-3 |
+# | :----- | :---- | -------: | -------: |
+# | Gurobi | false |       55 |       57 |
+# | Gurobi | true  |       12 |       17 |
+# | HiGHS  | false |      175 |      201 |
+# | HiGHS  | true  |      173 |      670 |
+# | HiGHS  | true* |       94 |      589 182, 201, 650 |
+
+# run_problem(
+#     optimizer_with_attributes(
+#         Gurobi.Optimizer,
+#         MOI.RelativeGapTolerance() => 1e-3,
+#     ),
+# )
+
+run_problem(
+    optimizer_with_attributes(
+        () -> MathOptLazy.Optimizer(Gurobi.Optimizer),
+        MOI.RelativeGapTolerance() => 1e-2,
+    ),
+)
+
+# run_problem(
+#     optimizer_with_attributes(
+#         HiGHS.Optimizer,
+#         MOI.RelativeGapTolerance() => 1e-3,
+#     ),
+# )
+
+# run_problem(
+#     optimizer_with_attributes(
+#         () -> MathOptLazy.Optimizer(HiGHS.Optimizer),
+#         MOI.RelativeGapTolerance() => 1e-3,
+#         "random_seed" => 123,
+#     ),
+# )
+
+# for (k, v) in model.internal.container.constraints
+#     if k isa POM.ConstraintKey{POM.FlowRateConstraint,PSY.Line}
+#         JuMP.set_attribute.(v, Gurobi.ConstraintAttribute("Lazy"), 1)
+#     end
+# end
+# @time solve = POM.solve!(model)
+
+# function solve_with_loop(model)
+#     jmp = POM.IOM.get_optimization_container(model).JuMPmodel
+#     constraints_lb, constraints_ub = Dict{Any,Any}(), Dict{Any,Any}()
+#     for (k, v) in model.internal.container.constraints
+#         if k == POM.ConstraintKey{POM.FlowRateConstraint,PSY.Line}("lb")
+#             for vi in v
+#                 constraints_lb[vi] = JuMP.constraint_object(vi)
+#             end
+#         elseif k == POM.ConstraintKey{POM.FlowRateConstraint,PSY.Line}("ub")
+#             for vi in v
+#                 constraints_ub[vi] = JuMP.constraint_object(vi)
+#             end
+#         end
+#     end
+#     JuMP.delete(jmp, [k for k in keys(constraints_lb)])
+#     JuMP.delete(jmp, [k for k in keys(constraints_ub)])
+#     JuMP.set_silent(jmp)
+#     total_solve_time = 0.0
+#     while true
+#         start_time = time()
+#         JuMP.optimize!(jmp)
+#         total_solve_time += time() - start_time
+#         n_constraints_added = 0
+#         for (k, c) in constraints_lb
+#             if JuMP.value(c.func) < c.set.lower
+#                 JuMP.@constraint(jmp, c.func in c.set)
+#                 n_constraints_added += 1
+#                 delete!(constraints_lb, k)
+#             end
+#         end
+#         for (k, c) in constraints_ub
+#             if JuMP.value(c.func) > c.set.upper
+#                 JuMP.@constraint(jmp, c.func in c.set)
+#                 n_constraints_added += 1
+#                 delete!(constraints_ub, k)
+#             end
+#         end
+#         if n_constraints_added == 0
+#             break
+#         else
+#             @show n_constraints_added
+#         end
+#     end
+#     @show total_solve_time
+#     return jmp
+# end
+
+# import PProf
+# POM.build!(model; output_dir = mktempdir(; cleanup = true))
+# PProf.@pprof POM.build!(model; output_dir = mktempdir(; cleanup = true))
+
+# using SnoopCompileCore
+# invs = @snoop_invalidations using PowerSystems, PowerOperationsModels
+# using SnoopCompile, AbstractTrees
+# trees = invalidation_trees(invs)