Constraints for multi-reservoir turbine outflow

[luke-kiernan]

If you attach 2 reservoirs to one HydroTurbine, it can draw up to outflow_limits.max from each reservoir, when I'd expect the limit to apply in aggregate. Is this a bug? Seems like it to me.

We even build the expression for total outflow for a hydro turbine here...but we don't add an upper bound of outflow_limits.max.

Minimum reproducible example, written by Claude

#=
Minimum reproducible example: HydroTurbine outflow_limits is applied per
(turbine, reservoir) edge rather than to the total turbine outflow.

When a turbine is connected to N head reservoirs, the solver can push up to
N * outflow_limits.max through it -- violating the turbine's aggregate capacity.

Run with:
  julia --project=test mre_turbine_outflow_bug.jl
=#

using PowerOperationsModels
using PowerSystems
using PowerSystemCaseBuilder
using InfrastructureSystems
using InfrastructureOptimizationModels
using InfrastructureSystems: SingleTimeSeries
using HiGHS
using DataFrames
using Dates

const PSY = PowerSystems
const PSB = PowerSystemCaseBuilder
const IS = InfrastructureSystems

sys = PSB.build_system(
    PSB.PSITestSystems,
    "c_sys5_hy_turbine_head";
    force_build = true,
    add_single_time_series = true,
)
transform_single_time_series!(sys, Dates.Hour(24), Dates.Hour(24))

turbine = only(get_components(HydroTurbine, sys))
res1    = only(get_components(HydroReservoir, sys))

# Clone the reservoir as a second head reservoir and wire it to the same turbine.
res2 = HydroReservoir(;
    name                  = "Water_Reservoir_2",
    available             = true,
    initial_level         = get_initial_level(res1),
    storage_level_limits  = get_storage_level_limits(res1),
    spillage_limits       = get_spillage_limits(res1),
    inflow                = get_inflow(res1),
    outflow               = get_outflow(res1),
    level_targets         = get_level_targets(res1),
    intake_elevation      = get_intake_elevation(res1),
    head_to_volume_factor = get_head_to_volume_factor(res1),
    level_data_type       = get_level_data_type(res1),
)
add_component!(sys, res2)
set_downstream_turbines!(res1, [turbine])
set_downstream_turbines!(res2, [turbine])

# Copy inflow time series from res1 onto res2.
inflow_ts = get_time_series(SingleTimeSeries, res1, "inflow")
add_time_series!(sys, res2, inflow_ts)
transform_single_time_series!(sys, Dates.Hour(24), Dates.Hour(24))

@assert length(get_connected_head_reservoirs(sys, turbine)) == 2

outflow_cap = PSY.get_outflow_limits(turbine).max
@info "Turbine outflow_limits.max (per-edge bound applied by code)" outflow_cap

template = OperationsProblemTemplate(NetworkModel(CopperPlatePowerModel))
set_device_model!(template, ThermalStandard, ThermalDispatchNoMin)
set_device_model!(template, PowerLoad, StaticPowerLoad)
set_device_model!(template, HydroTurbine, HydroTurbineWaterLinearDispatch)
set_device_model!(template,
    DeviceModel(HydroReservoir, HydroWaterModelReservoir))

model = DecisionModel(
    template, sys;
    name = "MRE",
    optimizer = optimizer_with_attributes(HiGHS.Optimizer),
    store_variable_names = true,
    optimizer_solve_log_print = false,
)
@assert build!(model; output_dir = mktempdir()) == ModelBuildStatus.BUILT
@assert solve!(model) == IS.Simulation.RunStatus.SUCCESSFULLY_FINALIZED

sol = OptimizationProblemOutputs(model)
total_flow_df = read_expression(sol, "TotalHydroFlowRateTurbineOutgoing__HydroTurbine")
total_flow = total_flow_df[!, "value"]

max_total = maximum(total_flow)
@info "Peak total outflow through turbine" max_total outflow_cap
if max_total > outflow_cap + 1e-6
    @warn "BUG CONFIRMED: turbine total outflow exceeds its outflow_limits.max" max_total outflow_cap
else
    @info "Bug NOT reproduced at these data values -- try increasing load / inflow."
end

result:

┌ Warning: BUG CONFIRMED: turbine total outflow exceeds its outflow_limits.max
│   max_total = 60.0
│   outflow_cap = 30.0
└ @ Main

[rodrigomha]

It could be. We have never tested this multi-reservoir into single turbine because is not common. Debugging it could be tricky.
Most of the existing hydro are based on one reservoir into multiple turbines, so we don't really have a test for this. I'm guessing the constraints for each flow are being created independently and we need to create a constraint that links both flows into a single one.