diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl
index fb941f4..97523ec 100644
--- a/src/PowerOperationsModels.jl
+++ b/src/PowerOperationsModels.jl
@@ -537,6 +537,7 @@ export HydroWaterFactorModel
 export HydroWaterModelReservoir
 export HydroTurbineBilinearDispatch
 export HydroTurbineWaterLinearDispatch
+export HydroTurbineMILPBilinearDispatch
 export HydroTurbineWaterLinearCommitment
 export HydroEnergyModelReservoir
 export HydroTurbineEnergyDispatch
diff --git a/src/core/formulations.jl b/src/core/formulations.jl
index 1e2d031..b353be1 100644
--- a/src/core/formulations.jl
+++ b/src/core/formulations.jl
@@ -326,6 +326,40 @@ Formulation type to add injection variables for a HydroTurbine connected to rese
 """
 struct HydroTurbineBilinearDispatch <: AbstractHydroDispatchFormulation end
 
+"""
+MILP formulation for the turbined-flow × head bilinear product in the hydro
+turbine power-output constraint. Adds injection variables for a HydroTurbine
+connected to reservoirs using a linearized approximation of the bilinear model.
+
+Selects the bilinear approximation scheme via `DeviceModel` attributes (so users
+do not need to depend on `InfrastructureOptimizationModels`). All attributes
+have defaults that reproduce the prior `Bin2` + `SolverSOS2`(4) behavior.
+
+# Attributes
+- `bilinear_approximation::String` (default `"bin2"`): top-level scheme.
+  Supported: `"bin2"`, `"hybs"`, `"nmdt"`, `"dnmdt"`, `"none"`.
+- `bilinear_quadratic_method::String` (default `"solver_sos2"`): inner quadratic
+  PWL method used when `bilinear_approximation ∈ {"bin2","hybs"}`. Supported:
+  `"solver_sos2"`, `"manual_sos2"`, `"sawtooth"`, `"epigraph"`, `"nmdt"`,
+  `"dnmdt"`, `"none"`.
+- `bilinear_n_segments::Int` (default `4`): PWL segment count or NMDT depth,
+  depending on the chosen method.
+- `bilinear_add_mccormick::Union{Bool,Nothing}` (default `nothing`): when
+  `nothing`, defers to the IOM struct default (`Bin2Config` → `true`,
+  `HybSConfig` → `false`). Ignored by `nmdt`/`dnmdt`/`none`.
+- `bilinear_epigraph_depth::Union{Int,Nothing}` (default `nothing`): when
+  `nothing`, defers to the IOM struct default (`NMDTQuadConfig` /
+  `DNMDTQuadConfig` → `3*depth`). `"hybs"` has no IOM default for this field
+  and *must* override.
+
+# Sentinel convention
+`nothing` means "use the IOM constructor's default value." POM does not
+duplicate IOM struct defaults; it just passes through the user's overrides.
+
+See: [`PowerSystems.HydroGen`](@extref).
+"""
+struct HydroTurbineMILPBilinearDispatch <: 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 +398,7 @@ These types share constructors.
 """
 const HydroTurbineWaterFormulation = Union{
     HydroTurbineBilinearDispatch,
+    HydroTurbineMILPBilinearDispatch,
     HydroTurbineWaterLinearDispatch,
     HydroTurbineWaterLinearCommitment,
 }
diff --git a/src/static_injector_models/hydro_generation.jl b/src/static_injector_models/hydro_generation.jl
index 93205d9..ec5a982 100644
--- a/src/static_injector_models/hydro_generation.jl
+++ b/src/static_injector_models/hydro_generation.jl
@@ -421,6 +421,41 @@ function get_default_attributes(
     return Dict{String, Any}("head_fraction_usage" => 0.0)
 end
 
+"""
+Default `DeviceModel` attributes for `HydroTurbineMILPBilinearDispatch`. The
+returned dictionary picks the bilinear approximation scheme used inside the
+turbine-power constraint; see [`HydroTurbineMILPBilinearDispatch`](@ref) for the
+full attribute reference and the `nothing` sentinel convention.
+
+The defaults reproduce the pre-rename behavior bit-for-bit:
+`IOM.Bin2Config(IOM.SolverSOS2QuadConfig(4))` with `add_mccormick = true`
+(inherited from `Bin2Config`'s one-argument constructor).
+"""
+function get_default_attributes(
+    ::Type{T},
+    ::Type{D},
+) where {T <: PSY.HydroTurbine, D <: HydroTurbineMILPBilinearDispatch}
+    return Dict{String, Any}(
+        # Top-level bilinear approximation scheme.
+        # Supported: "bin2", "hybs", "nmdt", "dnmdt", "none".
+        "bilinear_approximation" => "bin2",
+        # Inner quadratic PWL method (used when bilinear_approximation ∈ {"bin2","hybs"}).
+        # Supported: "solver_sos2", "manual_sos2", "sawtooth", "epigraph",
+        # "nmdt", "dnmdt", "none".
+        "bilinear_quadratic_method" => "solver_sos2",
+        # Segment count / discretization depth. SOS2/sawtooth/epigraph use it as
+        # n_segments; NMDT/DNMDT use it as depth.
+        "bilinear_n_segments" => 4,
+        # `nothing` ⇒ defer to IOM struct default
+        # (Bin2Config: true; HybSConfig: false; ignored otherwise).
+        "bilinear_add_mccormick" => nothing,
+        # `nothing` ⇒ defer to IOM struct default
+        # (HybSConfig has no default and must be overridden; NMDT/DNMDT
+        #  quadratic: 3*depth; ignored otherwise).
+        "bilinear_epigraph_depth" => nothing,
+    )
+end
+
 function get_default_attributes(
     ::Type{T},
     ::Type{D},
@@ -1817,6 +1852,177 @@ function add_constraints!(
     return
 end
 
+"""
+Build an `IOM.QuadraticApproxConfig` from attribute values.
+
+`n` is the segment count / discretization depth (`bilinear_n_segments`).
+`epi` is either an `Int` (`bilinear_epigraph_depth` override) or `nothing`,
+in which case IOM's struct default applies — relevant only for `NMDTQuadConfig`
+and `DNMDTQuadConfig`, which use `3*depth` when called with one argument.
+
+Errors with a list of supported method strings when `method` is unrecognized.
+"""
+function _build_quadratic_config(method::String, n::Int, epi)
+    if method == "solver_sos2"
+        return IOM.SolverSOS2QuadConfig(n)
+    elseif method == "manual_sos2"
+        return IOM.ManualSOS2QuadConfig(n)
+    elseif method == "sawtooth"
+        return IOM.SawtoothQuadConfig(n)
+    elseif method == "epigraph"
+        return IOM.EpigraphQuadConfig(n)
+    elseif method == "nmdt"
+        return epi === nothing ? IOM.NMDTQuadConfig(n) : IOM.NMDTQuadConfig(n, epi)
+    elseif method == "dnmdt"
+        return epi === nothing ? IOM.DNMDTQuadConfig(n) : IOM.DNMDTQuadConfig(n, epi)
+    elseif method == "none"
+        return IOM.NoQuadApproxConfig()
+    else
+        error(
+            "Unsupported bilinear_quadratic_method \"$(method)\". " *
+            "Supported: \"solver_sos2\", \"manual_sos2\", \"sawtooth\", " *
+            "\"epigraph\", \"nmdt\", \"dnmdt\", \"none\".",
+        )
+    end
+end
+
+"""
+Build an `IOM.BilinearApproxConfig` from the `DeviceModel`'s attributes.
+
+Reads `bilinear_approximation`, `bilinear_quadratic_method`,
+`bilinear_n_segments`, `bilinear_add_mccormick`, and `bilinear_epigraph_depth`
+from `model`. Sentinel (`nothing`) attribute values mean "let the IOM
+constructor's default apply" — POM never duplicates an IOM struct default.
+
+Errors when:
+- `bilinear_approximation` is not one of `"bin2"`, `"hybs"`, `"nmdt"`,
+  `"dnmdt"`, `"none"`.
+- `bilinear_approximation == "hybs"` but `bilinear_epigraph_depth === nothing`
+  (`HybSConfig` has no IOM-side default for `epigraph_depth`).
+
+See [`HydroTurbineMILPBilinearDispatch`](@ref) for the attribute reference.
+"""
+function _build_bilinear_config(
+    model::DeviceModel{<:PSY.HydroTurbine, HydroTurbineMILPBilinearDispatch},
+)
+    method = get_attribute(model, "bilinear_approximation")
+    n_segments = get_attribute(model, "bilinear_n_segments")
+    add_mc = get_attribute(model, "bilinear_add_mccormick")
+    epi = get_attribute(model, "bilinear_epigraph_depth")
+
+    if method == "bin2"
+        quad = _build_quadratic_config(
+            get_attribute(model, "bilinear_quadratic_method"),
+            n_segments,
+            epi,
+        )
+        return add_mc === nothing ? IOM.Bin2Config(quad) : IOM.Bin2Config(quad, add_mc)
+    elseif method == "hybs"
+        epi === nothing && error(
+            "bilinear_approximation = \"hybs\" requires a non-nothing " *
+            "bilinear_epigraph_depth attribute (IOM.HybSConfig has no default).",
+        )
+        quad = _build_quadratic_config(
+            get_attribute(model, "bilinear_quadratic_method"),
+            n_segments,
+            epi,
+        )
+        return if add_mc === nothing
+            IOM.HybSConfig(quad, epi)
+        else
+            IOM.HybSConfig(quad, epi, add_mc)
+        end
+    elseif method == "nmdt"
+        return IOM.NMDTBilinearConfig(n_segments)
+    elseif method == "dnmdt"
+        return IOM.DNMDTBilinearConfig(n_segments)
+    elseif method == "none"
+        return IOM.NoBilinearApproxConfig()
+    else
+        error(
+            "Unsupported bilinear_approximation \"$(method)\" for " *
+            "HydroTurbineMILPBilinearDispatch. " *
+            "Supported: \"bin2\", \"hybs\", \"nmdt\", \"dnmdt\", \"none\".",
+        )
+    end
+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 <: HydroTurbineMILPBilinearDispatch,
+    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!(
+            _build_bilinear_config(model),
+            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 ################################
 ############################################################################
diff --git a/src/static_injector_models/hydrogeneration_constructor.jl b/src/static_injector_models/hydrogeneration_constructor.jl
index 805c84a..6c3b2ea 100644
--- a/src/static_injector_models/hydrogeneration_constructor.jl
+++ b/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{HydroTurbineMILPBilinearDispatch},
+        Type{HydroTurbineWaterLinearDispatch},
+    },
 ) = nothing
 
 function construct_device!(
diff --git a/test/runtests.jl b/test/runtests.jl
index 83a9427..d290015 100644
--- a/test/runtests.jl
+++ b/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",
diff --git a/test/test_device_hydro_constructors.jl b/test/test_device_hydro_constructors.jl
index 2474198..f9b8653 100644
--- a/test/test_device_hydro_constructors.jl
+++ b/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 "HydroTurbineMILPBilinearDispatch: 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, HydroTurbineMILPBilinearDispatch)
+    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, HydroTurbineMILPBilinearDispatch)
+    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,