Bilinear delta formulation-based OptimizationContainer test

src/InfrastructureOptimizationModels.jl

@@ -481,7 +481,7 @@ export get_system_to_file, get_initialize_model, get_initialization_file
 export get_deserialize_initial_conditions, get_export_pwl_vars
 export get_check_numerical_bounds, get_allow_fails
 export get_optimizer_solve_log_print, get_calculate_conflict
-export get_detailed_optimizer_stats, get_direct_mode_optimizer
+export get_detailed_optimizer_stats, get_direct_model_optimizer
 export get_store_variable_names, get_export_optimization_model
 export use_time_series_cache
 export set_horizon!, set_initial_time!, set_warm_start!

src/core/optimization_container.jl

@@ -95,7 +95,7 @@ function OptimizationContainer(
         error("Default Time Series Type $T can't be abstract")
     end
 
-    if jump_model !== nothing && get_direct_mode_optimizer(settings)
+    if jump_model !== nothing && get_direct_model_optimizer(settings)
         throw(
             IS.ConflictingInputsError(
                 "Externally provided JuMP models are not compatible with the direct model keyword argument. Use JuMP.direct_model before passing the custom model",
@@ -302,7 +302,7 @@ function _finalize_jump_model!(container::OptimizationContainer, settings::Setti
         )
     end
 
-    if get_direct_mode_optimizer(settings)
+    if get_direct_model_optimizer(settings)
         optimizer = () -> MOI.instantiate(get_optimizer(settings))
         container.JuMPmodel = JuMP.direct_model(optimizer())
     elseif get_optimizer(settings) === nothing
@@ -332,18 +332,28 @@ function _finalize_jump_model!(container::OptimizationContainer, settings::Setti
     return
 end
 
+function intermediate_set_units_base_system!(sys::PSY.System, base)
+    PSY.set_units_base_system(sys, "SYSTEM_BASE")
+end
+
+function intermediate_get_forecast_initial_timestamp(sys::PSY.System)
+    return PSY.get_forecast_initial_timestamp(sys)
+end
+
 function init_optimization_container!(
     container::OptimizationContainer,
     network_model::NetworkModel{T},
-    sys::PSY.System,
+    sys::IS.InfrastructureSystemsContainer,
 ) where {T <: AbstractPowerModel}
-    PSY.set_units_base_system!(sys, "SYSTEM_BASE")
+    # PSY.set_units_base_system!(sys, "SYSTEM_BASE")
+    intermediate_set_units_base_system!(sys, "SYSTEM_BASE")
     # The order of operations matter
     settings = get_settings(container)
 
     if get_initial_time(settings) == UNSET_INI_TIME
         if get_default_time_series_type(container) <: PSY.AbstractDeterministic
-            set_initial_time!(settings, PSY.get_forecast_initial_timestamp(sys))
+            # set_initial_time!(settings, PSY.get_forecast_initial_timestamp(sys))
+            set_initial_time!(settings, intermediate_get_forecast_initial_timestamp(sys))
         elseif get_default_time_series_type(container) <: PSY.SingleTimeSeries
             ini_time, _ = PSY.check_time_series_consistency(sys, PSY.SingleTimeSeries)
             set_initial_time!(settings, ini_time)
@@ -437,7 +447,10 @@ end
 Execute the optimizer on the container's JuMP model, compute aux/dual variables,
 and return the run status. Called `solve_impl!(container, system)` in PSI.
 """
-function execute_optimizer!(container::OptimizationContainer, system::PSY.System)
+function execute_optimizer!(
+    container::OptimizationContainer,
+    system::IS.InfrastructureSystemsContainer,
+)
     optimizer_stats = get_optimizer_stats(container)
 
     jump_model = get_jump_model(container)

src/core/settings.jl

@@ -5,7 +5,7 @@ struct Settings
     warm_start::Base.RefValue{Bool}
     initial_time::Base.RefValue{Dates.DateTime}
     optimizer::Any  # Union{Nothing, MOI.OptimizerWithAttributes} or duck-typed optimizer
-    direct_mode_optimizer::Bool
+    direct_model_optimizer::Bool
     optimizer_solve_log_print::Bool
     detailed_optimizer_stats::Bool
     calculate_conflict::Bool
@@ -30,7 +30,7 @@ function Settings(
     horizon::Dates.Period = UNSET_HORIZON,
     resolution::Dates.Period = UNSET_RESOLUTION,
     optimizer = nothing,
-    direct_mode_optimizer::Bool = false,
+    direct_model_optimizer::Bool = false,
     optimizer_solve_log_print::Bool = false,
     detailed_optimizer_stats::Bool = false,
     calculate_conflict::Bool = false,
@@ -65,7 +65,7 @@ function Settings(
         Ref(warm_start),
         Ref(initial_time),
         optimizer_,
-        direct_mode_optimizer,
+        direct_model_optimizer,
         optimizer_solve_log_print,
         detailed_optimizer_stats,
         calculate_conflict,
@@ -148,7 +148,7 @@ get_allow_fails(settings::Settings) = settings.allow_fails
 get_optimizer_solve_log_print(settings::Settings) = settings.optimizer_solve_log_print
 get_calculate_conflict(settings::Settings) = settings.calculate_conflict
 get_detailed_optimizer_stats(settings::Settings) = settings.detailed_optimizer_stats
-get_direct_mode_optimizer(settings::Settings) = settings.direct_mode_optimizer
+get_direct_model_optimizer(settings::Settings) = settings.direct_model_optimizer
 get_store_variable_names(settings::Settings) = settings.store_variable_names
 get_rebuild_model(settings::Settings) = settings.rebuild_model
 get_export_optimization_model(settings::Settings) = settings.export_optimization_model

src/operation/decision_model.jl

@@ -49,7 +49,7 @@ Build the optimization problem of type M with the specific system and template.
   - `optimizer_solve_log_print::Bool = false`: Uses JuMP.unset_silent() to print the optimizer's log. By default all solvers are set to MOI.Silent()
   - `detailed_optimizer_stats::Bool = false`: True to save detailed optimizer stats log.
   - `calculate_conflict::Bool = false`: True to use solver to calculate conflicts for infeasible problems. Only specific solvers are able to calculate conflicts.
-  - `direct_mode_optimizer::Bool = false`: True to use the solver in direct mode. Creates a [JuMP.direct_model](https://jump.dev/JuMP.jl/dev/reference/models/#JuMP.direct_model).
+  - `direct_model_optimizer::Bool = false`: True to use the solver in direct mode. Creates a [JuMP.direct_model](https://jump.dev/JuMP.jl/dev/reference/models/#JuMP.direct_model).
   - `store_variable_names::Bool = false`: to store variable names in optimization model. Decreases the build times.
   - `rebuild_model::Bool = false`: It will force the rebuild of the underlying JuMP model with each call to update the model. It increases solution times, use only if the model can't be updated in memory.
   - `initial_time::Dates.DateTime = UNSET_INI_TIME`: Initial Time for the model solve.
@@ -112,7 +112,7 @@ function DecisionModel{M}(
     optimizer_solve_log_print = false,
     detailed_optimizer_stats = false,
     calculate_conflict = false,
-    direct_mode_optimizer = false,
+    direct_model_optimizer = false,
     store_variable_names = false,
     rebuild_model = false,
     export_optimization_model = false,
@@ -137,7 +137,7 @@ function DecisionModel{M}(
         calculate_conflict = calculate_conflict,
         optimizer_solve_log_print = optimizer_solve_log_print,
         detailed_optimizer_stats = detailed_optimizer_stats,
-        direct_mode_optimizer = direct_mode_optimizer,
+        direct_model_optimizer = direct_model_optimizer,
         check_numerical_bounds = check_numerical_bounds,
         store_variable_names = store_variable_names,
         rebuild_model = rebuild_model,

src/operation/emulation_model.jl

@@ -39,7 +39,7 @@ Build the optimization problem of type M with the specific system and template.
   - `calculate_conflict::Bool = false`: True to use solver to calculate conflicts for infeasible problems. Only specific solvers are able to calculate conflicts.
   - `optimizer_solve_log_print::Bool = false`: Uses JuMP.unset_silent() to print the optimizer's log. By default all solvers are set to MOI.Silent()
   - `detailed_optimizer_stats::Bool = false`: True to save detailed optimizer stats log.
-  - `direct_mode_optimizer::Bool = false`: True to use the solver in direct mode. Creates a [JuMP.direct_model](https://jump.dev/JuMP.jl/dev/reference/models/#JuMP.direct_model).
+  - `direct_model_optimizer::Bool = false`: True to use the solver in direct mode. Creates a [JuMP.direct_model](https://jump.dev/JuMP.jl/dev/reference/models/#JuMP.direct_model).
   - `store_variable_names::Bool = false`: True to store variable names in optimization model.
   - `rebuild_model::Bool = false`: It will force the rebuild of the underlying JuMP model with each call to update the model. It increases solution times, use only if the model can't be updated in memory.
   - `initial_time::Dates.DateTime = UNSET_INI_TIME`: Initial Time for the model solve.
@@ -106,7 +106,7 @@ function EmulationModel{M}(
     calculate_conflict = false,
     optimizer_solve_log_print = false,
     detailed_optimizer_stats = false,
-    direct_mode_optimizer = false,
+    direct_model_optimizer = false,
     check_numerical_bounds = true,
     store_variable_names = false,
     rebuild_model = false,
@@ -128,7 +128,7 @@ function EmulationModel{M}(
         calculate_conflict = calculate_conflict,
         optimizer_solve_log_print = optimizer_solve_log_print,
         detailed_optimizer_stats = detailed_optimizer_stats,
-        direct_mode_optimizer = direct_mode_optimizer,
+        direct_model_optimizer = direct_model_optimizer,
         check_numerical_bounds = check_numerical_bounds,
         store_variable_names = store_variable_names,
         rebuild_model = rebuild_model,

src/operation/operation_model_interface.jl

@@ -263,7 +263,7 @@ function _pre_solve_model_checks(model::OperationModel, optimizer = nothing)
             error("No Optimizer has been defined, can't solve the operational problem")
         end
     else
-        @assert get_direct_mode_optimizer(get_settings(model))
+        @assert get_direct_model_optimizer(get_settings(model))
     end
 
     optimizer_name = JuMP.solver_name(jump_model)

src/quadratic_approximations/nonlinear.jl

@@ -0,0 +1,27 @@
+struct BilinearProductExpression <: ExpressionType end
+
+function _add_bilinear!(
+    container::OptimizationContainer,
+    ::Type{C},
+    names::Vector{String},
+    time_steps::UnitRange{Int},
+    x_var_container,
+    y_var_container,
+    meta::String;
+) where {C <: IS.InfrastructureSystemsComponent}
+    z_container = add_expression_container!(
+        container,
+        BilinearProductExpression(),
+        C,
+        names,
+        time_steps;
+        expr_type = JuMP.QuadExpr,
+        meta,
+    )
+    for name in names, t in time_steps
+        z_expr = JuMP.QuadExpr()
+        JuMP.add_to_expression!(z_expr, x_var_container[name, t], y_var_container[name, t])
+        z_container[name, t] = z_expr
+    end
+    return z_container
+end

test/Project.toml

@@ -6,7 +6,9 @@ DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
+InfrastructureOptimizationModels = "bed98974-b02a-5e2f-9ee0-a103f5c45069"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -24,9 +26,9 @@ TimeSeries = "9e3dc215-6440-5c97-bce1-76c03772f85e"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
+[sources]
+InfrastructureSystems = {rev = "IS4", url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl"}
+
 [compat]
 HiGHS = "1"
 julia = "^1.10"
-
-[sources]
-InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"}

test/mocks/constructors.jl

@@ -3,6 +3,7 @@ Factory functions for quickly creating test fixtures.
 """
 
 using Dates
+using Random
 
 """
 Create a mock system with specified number of buses, generators, and loads.
@@ -72,7 +73,87 @@ function make_mock_thermal(
     bus = MockBus("bus1", 1, :PV),
     limits = (min = 0.0, max = 100.0),
     base_power = 100.0,
-    operation_cost = MockOperationCost(0.0),
+    operation_cost = MockProportionalCost(0.0),
 )
     return MockThermalGen(name, available, bus, limits, base_power, operation_cost)
 end
+
+"""
+Generate convex piecewise-linear cost curve points with `n_tranches` segments over [0, pmax].
+Returns a vector of (x, y) tuples with strictly increasing slopes.
+"""
+function _random_convex_pwl_points(n_tranches::Int, pmax::Float64, rng)
+    xs = sort(rand(rng, n_tranches - 1)) .* pmax
+    points = [(0.0, 0.0)]
+    cumulative_cost = 0.0
+    prev_x = 0.0
+    slope = 5.0 + 20.0 * rand(rng)
+    for x in xs
+        cumulative_cost += slope * (x - prev_x)
+        push!(points, (x, cumulative_cost))
+        prev_x = x
+        slope += 5.0 + 10.0 * rand(rng)
+    end
+    cumulative_cost += slope * (pmax - prev_x)
+    push!(points, (pmax, cumulative_cost))
+    return points
+end
+
+"""
+Create a mock system on `n` nodes where the graph is connected (and fairly random),
+half the nodes have generators, and half have loads.
+"""
+function make_mock_test_network(n::Int; max_tranches::Int = 4, seed::Int = 42)
+    @assert n >= 2 "Need at least 2 nodes for a network"
+    rng = MersenneTwister(seed)
+    sys = MockSystem(1.0)
+
+    # Create buses
+    buses = [MockBus("bus$i", i, :PV) for i in 1:n]
+    for bus in buses
+        add_component!(sys, bus)
+    end
+
+    # Build a connected graph: chain 1-2-3-...-n plus random extra edges
+    branch_pairs = Set{Tuple{Int, Int}}()
+    perm = shuffle(rng, 1:n)
+    for i in 1:(n - 1)
+        push!(branch_pairs, (perm[i], perm[i + 1]))
+    end
+    # Add some random cross-links for variety
+    n_extra = div(n, 3)
+    for _ in 1:n_extra
+        a, b = rand(rng, 1:n), rand(rng, 1:n)
+        if a != b
+            pair = a < b ? (a, b) : (b, a)
+            push!(branch_pairs, pair)
+        end
+    end
+    for (idx, (i, j)) in enumerate(branch_pairs)
+        r = 0.005 + 0.005 * rand(rng)
+        branch = MockBranch("branch$idx", true, buses[i], buses[j], 1.0, r)
+        add_component!(sys, branch)
+    end
+
+    # Half the nodes get generators with convex PWL costs, the other half get loads
+    n_gens = div(n, 2)
+    for i in 1:n_gens
+        n_tranches = rand(rng, 2:max_tranches)
+        pmax = 1.5 * rand(rng)
+        points = _random_convex_pwl_points(n_tranches, pmax, rng)
+
+        pwl = IS.PiecewiseLinearData(points)
+        cost_curve = IS.CostCurve(IS.InputOutputCurve(pwl))
+        op_cost = MockOperationalCost(cost_curve, 0.0, 0.0)
+        gen = MockThermalGen(
+            "gen$i", true, buses[i], (min = 0.0, max = pmax), 1.0, op_cost,
+        )
+        add_component!(sys, gen)
+    end
+    for i in (n_gens + 1):n
+        load = MockLoad("load$(i - n_gens)", true, buses[i], 0.5)
+        add_component!(sys, load)
+    end
+
+    return sys
+end