Add hydro formulation with bilinear approximation

src/PowerOperationsModels.jl

@@ -537,6 +537,7 @@ export HydroWaterFactorModel
 export HydroWaterModelReservoir
 export HydroTurbineBilinearDispatch
 export HydroTurbineWaterLinearDispatch
+export HydroTurbineBin2BilinearDispatch
 export HydroTurbineWaterLinearCommitment
 export HydroEnergyModelReservoir
 export HydroTurbineEnergyDispatch

src/core/formulations.jl

@@ -326,6 +326,11 @@ Formulation type to add injection variables for a HydroTurbine connected to rese
 """
 struct HydroTurbineBilinearDispatch <: AbstractHydroDispatchFormulation end
 
+"""
+Formulation type to add injection variables for a HydroTurbine connected to reservoirs using a bilinear model (with water flow variables) [`PowerSystems.HydroGen`](@extref). Uses a linearized approximation.
+"""
+struct HydroTurbineBin2BilinearDispatch <: AbstractHydroDispatchFormulation end
+
 """
 Formulation type to add injection variables for a HydroTurbine connected to reservoirs using a linear model [`PowerSystems.HydroGen`](@extref).
 The model assumes a shallow reservoir. The head for the conversion between water flow and power can be approximated as a linear function of the water flow on which the head elevation is always the intake elevation.
@@ -364,6 +369,7 @@ These types share constructors.
 """
 const HydroTurbineWaterFormulation = Union{
     HydroTurbineBilinearDispatch,
+    HydroTurbineBin2BilinearDispatch,
     HydroTurbineWaterLinearDispatch,
     HydroTurbineWaterLinearCommitment,
 }

src/static_injector_models/hydro_generation.jl

@@ -1817,6 +1817,82 @@ function add_constraints!(
     return
 end
 
+"""
+This function define the relationship between turbined flow and power produced with a linear approximation for the bilinear product.
+"""
+function add_constraints!(
+    container::OptimizationContainer,
+    sys::PSY.System,
+    ::Type{TurbinePowerOutputConstraint},
+    devices::IS.FlattenIteratorWrapper{V},
+    model::DeviceModel{V, W},
+    ::NetworkModel{X},
+) where {
+    V <: PSY.HydroTurbine,
+    W <: HydroTurbineBin2BilinearDispatch,
+    X <: PM.AbstractPowerModel,
+}
+    time_steps = get_time_steps(container)
+    base_power = get_model_base_power(container)
+    names = PSY.get_name.(devices)
+    constraint =
+        add_constraints_container!(
+            container,
+            TurbinePowerOutputConstraint,
+            V,
+            names,
+            time_steps,
+        )
+    power = get_variable(container, ActivePowerVariable, V)
+    flow = get_variable(container, HydroTurbineFlowRateVariable, V)
+    head = get_variable(container, HydroReservoirHeadVariable, PSY.HydroReservoir)
+    for d in devices
+        name = PSY.get_name(d)
+        conversion_factor = PSY.get_conversion_factor(d)
+        reservoirs = filter(PSY.get_available, PSY.get_connected_head_reservoirs(sys, d))
+        powerhouse_elevation = PSY.get_powerhouse_elevation(d)
+
+        fh_prod = IOM._add_bilinear_approx!(
+            IOM.Bin2Config(IOM.SolverSOS2QuadConfig(4)),
+            container,
+            V,
+            PSY.get_name.(reservoirs),
+            time_steps,
+            flow[name, :, :],
+            head,
+            repeat(
+                [
+                    (
+                    min = get_variable_lower_bound(HydroTurbineFlowRateVariable, d, W),
+                    max = get_variable_upper_bound(HydroTurbineFlowRateVariable, d, W),
+                )
+                ], length(reservoirs)),
+            [
+                (
+                    min = get_variable_lower_bound(HydroReservoirHeadVariable, res, W),
+                    max = get_variable_upper_bound(HydroReservoirHeadVariable, res, W),
+                ) for res in reservoirs
+            ],
+            "$(get_name(d))_FlowHeadProduct",
+        )
+
+        for t in time_steps
+            constraint[name, t] = JuMP.@constraint(
+                container.JuMPmodel,
+                power[name, t] ==
+                GRAVITATIONAL_CONSTANT * WATER_DENSITY * conversion_factor *
+                sum(
+                    fh_prod[PSY.get_name(res), t]
+                    -
+                    powerhouse_elevation * flow[name, PSY.get_name(res), t]
+                    for res in reservoirs
+                ) / (1e6 * base_power)
+            )
+        end
+    end
+    return
+end
+
 ############################################################################
 ############################### Expressions ################################
 ############################################################################

src/static_injector_models/hydrogeneration_constructor.jl

@@ -1809,7 +1809,11 @@ _maybe_add_on_variables!(
 _maybe_add_on_variables!(
     ::OptimizationContainer,
     devices,
-    ::Union{Type{HydroTurbineBilinearDispatch}, Type{HydroTurbineWaterLinearDispatch}},
+    ::Union{
+        Type{HydroTurbineBilinearDispatch},
+        Type{HydroTurbineBin2BilinearDispatch},
+        Type{HydroTurbineWaterLinearDispatch},
+    },
 ) = nothing
 
 function construct_device!(

test/runtests.jl

@@ -22,11 +22,15 @@ const DISABLED_TEST_FILES = [  # Can generate with ls -1 test | grep "test_.*.jl
 # "test_device_synchronous_condenser_constructors.jl",
 # "test_device_thermal_generation_constructors.jl",
 # "test_formulation_combinations.jl",
+# "test_import_export_cost.jl",
 # "test_initialization_problem.jl",
+# "test_is_time_variant_proportional.jl",
+# "test_market_bid_cost.jl",
+# "test_mbc_parameter_population.jl",
 # "test_model_decision.jl",
 # "test_multi_interval.jl",
-# "test_network_constructors_with_dlr.jl",
 # "test_network_constructors.jl",
+# "test_network_constructors_with_dlr.jl",
 # "test_problem_template.jl",
 # "test_storage_device_models.jl",
 # "test_transfer_initial_conditions.jl",

test/test_device_hydro_constructors.jl

@@ -678,11 +678,204 @@ end
     total_outflow = sum(df_outflow[!, :value])
     total_spillage = sum(hydro_spillage_df[!, :value])
 
+    # Tolerance covers accumulated rounding in the m^3/s → km^3 unit conversion
+    # plus the MILP bilinear approximation; water-balance closure within 1e-4 km^3
+    # is well inside HiGHS' default feasibility tolerance for this problem size.
+    tol = 1e-4
     calculated_vf =
         (hydro_vol_df[1, :value]) +
-        ((total_inflow - total_outflow - total_spillage) * 3600 * 1e-9)
+        (
+            (total_inflow - total_outflow - total_spillage) *
+            POM.SECONDS_IN_HOUR *
+            POM.M3_TO_KM3
+        )
+
+    @test isapprox(calculated_vf, hydro_vol_df[end, :value], atol = tol)
+
+    psi_checksolve_test(
+        model,
+        [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL, MOI.LOCALLY_SOLVED],
+        210949.49,
+        1000,
+    )
+end
+
+@testset "HydroTurbineBin2BilinearDispatch: variable-bound plumbing to IOM" begin
+    # Spot-check that POM forwards PSY device data to JuMP without unit conversion.
+    # Outflow limits are m^3/s and storage_level_limits is meters (HEAD reservoir),
+    # so JuMP variables should carry those values verbatim.
+    output_dir = mktempdir(; cleanup = true)
+
+    sys = PSB.build_system(PSITestSystems, "c_sys5_hy_turbine_head")
+    template = OperationsProblemTemplate()
+    set_device_model!(template, HydroTurbine, HydroTurbineBin2BilinearDispatch)
+    set_device_model!(template, HydroReservoir, HydroWaterModelReservoir)
+    set_device_model!(template, ThermalStandard, ThermalDispatchNoMin)
+    set_device_model!(template, PowerLoad, StaticPowerLoad)
+
+    model = DecisionModel(
+        template,
+        sys;
+        optimizer = HiGHS_optimizer,
+        store_variable_names = true,
+    )
+    @test build!(model; output_dir = output_dir) == ModelBuildStatus.BUILT
+
+    container = IOM.get_optimization_container(model)
+    time_steps = IOM.get_time_steps(container)
+    flow = IOM.get_variable(container, HydroTurbineFlowRateVariable, HydroTurbine)
+    head = IOM.get_variable(container, HydroReservoirHeadVariable, HydroReservoir)
+
+    reservoirs = collect(get_components(HydroReservoir, sys))
+    @test !isempty(reservoirs)
+    for res in reservoirs
+        # HEAD reservoirs store level limits directly in meters; bounds should match.
+        @test PSY.get_level_data_type(res) == PSY.ReservoirDataType.HEAD
+        limits = PSY.get_storage_level_limits(res)
+        r_name = PSY.get_name(res)
+        for t in time_steps
+            v = head[r_name, t]
+            @test JuMP.lower_bound(v) == limits.min
+            @test JuMP.upper_bound(v) == limits.max
+        end
+    end
+
+    turbines = collect(get_components(HydroTurbine, sys))
+    @test !isempty(turbines)
+    for turbine in turbines
+        outflow = PSY.get_outflow_limits(turbine)
+        @test outflow !== nothing
+        t_name = PSY.get_name(turbine)
+        for res in reservoirs, t in time_steps
+            v = flow[t_name, PSY.get_name(res), t]
+            @test JuMP.lower_bound(v) == outflow.min
+            @test JuMP.upper_bound(v) == outflow.max
+        end
+    end
+end
+
+@testset "HydroTurbineBilinearDispatch: TurbinePowerOutputConstraint unit conversion" begin
+    # Spot-check that POM's per-unit conversion (g*ρ*conv_factor / (1e6 * base_power))
+    # lands in JuMP exactly as expected. The pure bilinear formulation produces a clean
+    # quadratic constraint whose coefficients are easy to read off.
+    output_dir = mktempdir(; cleanup = true)
+
+    sys = PSB.build_system(PSITestSystems, "c_sys5_hy_turbine_head")
+    template = OperationsProblemTemplate()
+    set_device_model!(template, HydroTurbine, HydroTurbineBilinearDispatch)
+    set_device_model!(template, HydroReservoir, HydroWaterModelReservoir)
+    set_device_model!(template, ThermalStandard, ThermalDispatchNoMin)
+    set_device_model!(template, PowerLoad, StaticPowerLoad)
+
+    model = DecisionModel(
+        template,
+        sys;
+        optimizer = ipopt_optimizer,
+        store_variable_names = true,
+    )
+    @test build!(model; output_dir = output_dir) == ModelBuildStatus.BUILT
+
+    container = IOM.get_optimization_container(model)
+    time_steps = IOM.get_time_steps(container)
+    base_power = IOM.get_model_base_power(container)
+    constraint = IOM.get_constraint(container, TurbinePowerOutputConstraint, HydroTurbine)
+    power = IOM.get_variable(container, ActivePowerVariable, HydroTurbine)
+    flow = IOM.get_variable(container, HydroTurbineFlowRateVariable, HydroTurbine)
+    head = IOM.get_variable(container, HydroReservoirHeadVariable, HydroReservoir)
+
+    for turbine in get_components(HydroTurbine, sys)
+        t_name = PSY.get_name(turbine)
+        conv = PSY.get_conversion_factor(turbine)
+        elev = PSY.get_powerhouse_elevation(turbine)
+        reservoirs =
+            filter(PSY.get_available, PSY.get_connected_head_reservoirs(sys, turbine))
+        scale = POM.GRAVITATIONAL_CONSTANT * POM.WATER_DENSITY * conv / (1e6 * base_power)
+
+        # power = scale * Σ_res (head_res * flow - elev * flow)
+        #   ⇒ rearranged: power - scale * Σ head*flow + scale*elev * Σ flow == 0
+        for t in time_steps
+            con = constraint[t_name, t]
+            con_obj = JuMP.constraint_object(con)
+            @test con_obj.set isa MOI.EqualTo{Float64}
+
+            # Linear coefficients: power has +1, each flow has +scale*elev.
+            @test JuMP.normalized_coefficient(con, power[t_name, t]) ≈ 1.0
+            for res in reservoirs
+                @test JuMP.normalized_coefficient(
+                    con,
+                    flow[t_name, PSY.get_name(res), t],
+                ) ≈ scale * elev
+
+                # Quadratic cross term head*flow has -scale.
+                quad_terms = con_obj.func.terms
+                pair = JuMP.UnorderedPair(
+                    head[PSY.get_name(res), t],
+                    flow[t_name, PSY.get_name(res), t],
+                )
+                @test get(quad_terms, pair, 0.0) ≈ -scale
+            end
+        end
+    end
+end
+
+@testset "Solve HydroWaterModelReservoir with bilinear approximations" begin
+    output_dir = mktempdir(; cleanup = true)
+
+    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy_turbine_head")
+    reservoir = only(get_components(HydroReservoir, c_sys5_hy))
+    hydro_inflow_ts = get_time_series_array(Deterministic, reservoir, "inflow")
+
+    template = OperationsProblemTemplate()
+    set_device_model!(template, HydroTurbine, HydroTurbineBin2BilinearDispatch)
+    set_device_model!(template, HydroReservoir, HydroWaterModelReservoir)
+
+    set_device_model!(template, ThermalStandard, ThermalDispatchNoMin)
+    set_device_model!(template, PowerLoad, StaticPowerLoad)
+
+    model = DecisionModel(
+        template,
+        c_sys5_hy;
+        optimizer = HiGHS_optimizer,
+        store_variable_names = true,
+    )
+
+    @test build!(model; output_dir = output_dir) ==
+          ModelBuildStatus.BUILT
+
+    @test solve!(model; output_dir = output_dir) ==
+          IS.Simulation.RunStatus.SUCCESSFULLY_FINALIZED
+
+    outputs = OptimizationProblemOutputs(model)
+
+    psi_checkobjfun_test(model, AffExpr)
+
+    df_outflow = read_expression(outputs, "TotalHydroFlowRateTurbineOutgoing__HydroTurbine")
+    hydro_vol_df =
+        read_variables(outputs, [(HydroReservoirVolumeVariable, HydroReservoir)])["HydroReservoirVolumeVariable__HydroReservoir"]
+    hydro_head_df =
+        read_variables(outputs, [(HydroReservoirHeadVariable, HydroReservoir)])["HydroReservoirHeadVariable__HydroReservoir"]
+    hydro_spillage_df =
+        read_variables(outputs, [(WaterSpillageVariable, HydroReservoir)])["WaterSpillageVariable__HydroReservoir"]
+    hydro_inflow_df =
+        read_parameters(outputs, [(InflowTimeSeriesParameter, HydroReservoir)])["InflowTimeSeriesParameter__HydroReservoir"]
+
+    total_inflow = sum(values(hydro_inflow_ts))
+    total_outflow = sum(df_outflow[!, :value])
+    total_spillage = sum(hydro_spillage_df[!, :value])
+
+    # Tolerance covers accumulated rounding in the m^3/s → km^3 unit conversion
+    # plus the MILP bilinear approximation; water-balance closure within 1e-4 km^3
+    # is well inside HiGHS' default feasibility tolerance for this problem size.
+    tol = 1e-4
+    calculated_vf =
+        (hydro_vol_df[1, :value]) +
+        (
+            (total_inflow - total_outflow - total_spillage) *
+            POM.SECONDS_IN_HOUR *
+            POM.M3_TO_KM3
+        )
 
-    @test abs(calculated_vf - hydro_vol_df[end, :value]) <= 1e-4
+    @test isapprox(calculated_vf, hydro_vol_df[end, :value], atol = tol)
 
     psi_checksolve_test(
         model,