HydroTurbine + Reservoir for WaterModel

Note

HydroPowerSimulations.jl is an extension library of PowerSimulations.jl for modeling hydro units. Users are encouraged to review the single-step tutorial in PowerSimulations.jl before this tutorial.

Load packages

julia> using PowerSystems
julia> using PowerSimulations
julia> using HydroPowerSimulations
julia> using PowerSystemCaseBuilder
julia> using Ipopt # solver

Data

Note

PowerSystemCaseBuilder.jl is a helper library that makes it easier to reproduce examples in the documentation and tutorials. Normally you would pass your local files to create the system data instead of calling the function build_system. For more details visit PowerSystemCaseBuilder README

julia> sys = build_system(PSITestSystems, "c_sys5_hy_turbine_head")┌ Info: Building new system c_sys5_hy_turbine_head from raw data
  sys_descriptor.raw_data = "/home/runner/.julia/artifacts/edcb5940e84a802a86ad4f2223214d33121ac044/PowerSystemsTestData-4.0.2/psy_data/data_5bus_pu.jl"
[ Info: Serialized time series data to /home/runner/.julia/packages/PowerSystemCaseBuilder/zW01F/data/serialized_system/4e67b70ea6977dbe21c7731d72cdc1494adf072a7f3f08d921db740cf264ce79/c_sys5_hy_turbine_head_time_series_storage.h5.
[ Info: Serialized System to /home/runner/.julia/packages/PowerSystemCaseBuilder/zW01F/data/serialized_system/4e67b70ea6977dbe21c7731d72cdc1494adf072a7f3f08d921db740cf264ce79/c_sys5_hy_turbine_head.json
[ Info: Serialized System metadata to /home/runner/.julia/packages/PowerSystemCaseBuilder/zW01F/data/serialized_system/4e67b70ea6977dbe21c7731d72cdc1494adf072a7f3f08d921db740cf264ce79/c_sys5_hy_turbine_head_metadata.json
              System
┌───────────────────┬─────────────┐
│ Property           Value       │
├───────────────────┼─────────────┤
│ Name              │             │
│ Description       │             │
│ System Units Base │ SYSTEM_BASE │
│ Base Power        │ 100.0       │
│ Base Frequency    │ 60.0        │
│ Num Components    │ 27          │
└───────────────────┴─────────────┘

     Static Components
┌─────────────────┬───────┐
│ Type             Count │
├─────────────────┼───────┤
│ ACBus           │ 5     │
│ Arc             │ 6     │
│ HydroReservoir  │ 1     │
│ HydroTurbine    │ 1     │
│ Line            │ 6     │
│ PowerLoad       │ 3     │
│ ThermalStandard │ 5     │
└─────────────────┴───────┘

                                Forecast Summary
┌────────────────┬────────────────┬──────────────────┬──────────────────┬───────
│ owner_type      owner_category  name              time_series_type  init ⋯
│ String          String          String            String            Stri ⋯
├────────────────┼────────────────┼──────────────────┼──────────────────┼───────
│ HydroReservoir │ Component      │ inflow           │ Deterministic    │ 2024 ⋯
│ HydroReservoir │ Component      │ outflow          │ Deterministic    │ 2024 ⋯
│ PowerLoad      │ Component      │ max_active_power │ Deterministic    │ 2024 ⋯
└────────────────┴────────────────┴──────────────────┴──────────────────┴───────
                                                               6 columns omitted

With a single PowerSystems.HydroTurbine connected downstream to a PowerSystems.HydroReservoir:

julia> hy = only(get_components(HydroTurbine, sys))HydroTurbine: Water_Turbine:
   name: Water_Turbine
   available: true
   bus: ACBus: nodeC
   active_power: 0.0
   reactive_power: 0.0
   rating: 5.2
   active_power_limits: (min = 0.0, max = 5.2)
   reactive_power_limits: (min = -3.9, max = 3.9)
   base_power: 100.0
   operation_cost: PowerSystems.HydroGenerationCost composed of variable: InfrastructureSystems.CostCurve{InfrastructureSystems.LinearCurve}
   powerhouse_elevation: 317.12
   ramp_limits: nothing
   time_limits: nothing
   outflow_limits: (min = 0.0, max = 30.0)
   efficiency: 1.0
   turbine_type: PowerSystems.HydroTurbineTypeModule.HydroTurbineType.UNKNOWN = 0
   conversion_factor: 1.0
   prime_mover_type: PowerSystems.PrimeMoversModule.PrimeMovers.HY = 16
   travel_time: nothing
   services: 0-element Vector{PowerSystems.Service}
   dynamic_injector: nothing
   ext: Dict{String, Any}()
   InfrastructureSystems.SystemUnitsSettings:
      base_value: 100.0
      unit_system: InfrastructureSystems.UnitSystemModule.UnitSystem.SYSTEM_BASE = 0
   has_supplemental_attributes: false
   has_time_series: false
julia> res = only(get_components(HydroReservoir, sys))HydroReservoir: Water_Reservoir:
   name: Water_Reservoir
   available: true
   storage_level_limits: (min = 463.5, max = 555.5)
   initial_level: 0.9
   spillage_limits: nothing
   inflow: 1.0
   outflow: 0.0
   level_targets: 1.0
   intake_elevation: 463.3
   head_to_volume_factor: InfrastructureSystems.LinearCurve(302376.2, 0.0)
   upstream_turbines: 0-element Vector{PowerSystems.HydroUnit}
   downstream_turbines: 1-element Vector{PowerSystems.HydroUnit}
   upstream_reservoirs: 0-element Vector{PowerSystems.Device}
   operation_cost:
   level_data_type: PowerSystems.ReservoirDataTypeModule.ReservoirDataType.HEAD = 3
   ext: Dict{String, Any}()
   InfrastructureSystems.SystemUnitsSettings:
      base_value: 100.0
      unit_system: InfrastructureSystems.UnitSystemModule.UnitSystem.SYSTEM_BASE = 0
   has_supplemental_attributes: false
   has_time_series: true

Note that the reservoir has a level_data_type of HEAD, that implies its storage level limits data are in meters (above the sea level) and refer to the hydraulic head levels. That means that its available capacity lies with its hydraulic head being within 463.5 and 555.5 meters, and its intake elevation is at 463.3 meters. In addition note that the elevation of the turbine is on 317.12 meters above the sea level.

Decision Model

Setting up the formulations based on PowerSimulations.jl:

julia> template = ProblemTemplate(PTDFPowerModel)                     Network Model
┌────────────────────┬─────────────────────────────────┐
│ Network Model      │ PowerSimulations.PTDFPowerModel │
│ Slacks             │ false                           │
│ PTDF               │ false                           │
│ Duals              │ None                            │
│ HVDC Network Model │ None                            │
└────────────────────┴─────────────────────────────────┘

            Device Models
┌─────────────┬─────────────┬────────┐
│ Device Type  Formulation  Slacks │
└─────────────┴─────────────┴────────┘
julia> set_device_model!(template, ThermalStandard, ThermalBasicDispatch)
julia> set_device_model!(template, PowerLoad, StaticPowerLoad)
julia> set_device_model!(template, Line, StaticBranch)

but, now we also include the HydroTurbine using HydroTurbineBilinearDispatch:

julia> set_device_model!(template, HydroTurbine, HydroTurbineBilinearDispatch)

This is a nonlinear model that to compute its output power requires the bilinear term head times water flow. For that purpose the non-convex Ipopt solver will be used to solve this problem.

In addition, we need to use the water model for the HydroReservoir via HydroWaterModelReservoir.

julia> set_device_model!(template, HydroReservoir, HydroWaterModelReservoir)

With the template properly set-up, we construct, build and solve the optimization problem:

julia> model = DecisionModel(template, sys; optimizer = Ipopt.Optimizer)[ Info: Overriding time_series_cache_size because time series is stored in memory
                     Network Model
┌────────────────────┬─────────────────────────────────┐
│ Network Model      │ PowerSimulations.PTDFPowerModel │
│ Slacks             │ false                           │
│ PTDF               │ false                           │
│ Duals              │ None                            │
│ HVDC Network Model │ None                            │
└────────────────────┴─────────────────────────────────┘

                                 Device Models
┌──────────────────────────────┬───────────────────────────────────────┬────────
│ Device Type                   Formulation                            Slack ⋯
├──────────────────────────────┼───────────────────────────────────────┼────────
│ PowerSystems.ThermalStandard │ PowerSimulations.ThermalBasicDispatch │ false ⋯
│ PowerSystems.HydroReservoir  │ HydroWaterModelReservoir              │ false ⋯
│ PowerSystems.PowerLoad       │ PowerSimulations.StaticPowerLoad      │ false ⋯
│ PowerSystems.HydroTurbine    │ HydroTurbineBilinearDispatch          │ false ⋯
└──────────────────────────────┴───────────────────────────────────────┴────────
                                                                1 column omitted

                        Branch Models
┌───────────────────┬───────────────────────────────┬────────┐
│ Branch Type        Formulation                    Slacks │
├───────────────────┼───────────────────────────────┼────────┤
│ PowerSystems.Line │ PowerSimulations.StaticBranch │ false  │
└───────────────────┴───────────────────────────────┴────────┘
julia> build!(model; output_dir = mktempdir())InfrastructureSystems.Optimization.ModelBuildStatusModule.ModelBuildStatus.BUILT = 0
julia> solve!(model) ****************************************************************************** This program contains Ipopt, a library for large-scale nonlinear optimization. Ipopt is released as open source code under the Eclipse Public License (EPL). For more information visit https://github.com/coin-or/Ipopt ****************************************************************************** InfrastructureSystems.Simulation.RunStatusModule.RunStatus.SUCCESSFULLY_FINALIZED = 0

Exploring Results

Results can be explored using:

julia> res = OptimizationProblemResults(model)Start: 2024-01-01T00:00:00
End: 2024-01-01T23:00:00
Resolution: 60 minutes

     PowerSimulations Problem Expressions Results
┌─────────────────────────────────────────────────────┐
│ ProductionCostExpression__HydroTurbine              │
│ PTDFBranchFlow__Line                                │
│ ActivePowerBalance__ACBus                           │
│ ProductionCostExpression__ThermalStandard           │
│ TotalHydroFlowRateReservoirIncoming__HydroReservoir │
│ ActivePowerBalance__System                          │
│ TotalHydroFlowRateReservoirOutgoing__HydroReservoir │
│ TotalHydroFlowRateTurbineOutgoing__HydroTurbine     │
└─────────────────────────────────────────────────────┘

 PowerSimulations Problem Parameters Results
┌────────────────────────────────────────────┐
│ ActivePowerTimeSeriesParameter__PowerLoad  │
│ InflowTimeSeriesParameter__HydroReservoir  │
│ OutflowTimeSeriesParameter__HydroReservoir │
└────────────────────────────────────────────┘

   PowerSimulations Problem Variables Results
┌──────────────────────────────────────────────┐
│ ActivePowerVariable__ThermalStandard         │
│ ActivePowerVariable__HydroTurbine            │
│ HydroReservoirHeadVariable__HydroReservoir   │
│ HydroWaterSurplusVariable__HydroReservoir    │
│ HydroWaterShortageVariable__HydroReservoir   │
│ HydroReservoirVolumeVariable__HydroReservoir │
│ WaterSpillageVariable__HydroReservoir        │
└──────────────────────────────────────────────┘

Use read_variable to read in the dispatch variable results for the hydro:

julia> var =
           read_variable(res, "ActivePowerVariable__HydroTurbine"; table_format = TableFormat.WIDE)24×2 DataFrame
 Row  DateTime             Water_Turbine 
      DateTime             Float64?      
─────┼────────────────────────────────────
   1 │ 2024-01-01T00:00:00        53.7386
   2 │ 2024-01-01T01:00:00        53.6704
   3 │ 2024-01-01T02:00:00        53.6022
   4 │ 2024-01-01T03:00:00        53.534
   5 │ 2024-01-01T04:00:00        53.4658
   6 │ 2024-01-01T05:00:00        53.3976
   7 │ 2024-01-01T06:00:00        53.3293
   8 │ 2024-01-01T07:00:00        53.2611
  ⋮  │          ⋮                 ⋮
  18 │ 2024-01-01T17:00:00        52.5789
  19 │ 2024-01-01T18:00:00        52.5107
  20 │ 2024-01-01T19:00:00        52.4425
  21 │ 2024-01-01T20:00:00        52.3743
  22 │ 2024-01-01T21:00:00        52.306
  23 │ 2024-01-01T22:00:00        52.2378
  24 │ 2024-01-01T23:00:00        52.1696
                            9 rows omitted

or the water flowing through the turbine (in m³/s):

julia> var = read_expression(
           res,
           "TotalHydroFlowRateTurbineOutgoing__HydroTurbine";
           table_format = TableFormat.WIDE,
       )24×2 DataFrame
 Row  DateTime             Water_Turbine 
      DateTime             Float64?      
─────┼────────────────────────────────────
   1 │ 2024-01-01T00:00:00           30.0
   2 │ 2024-01-01T01:00:00           30.0
   3 │ 2024-01-01T02:00:00           30.0
   4 │ 2024-01-01T03:00:00           30.0
   5 │ 2024-01-01T04:00:00           30.0
   6 │ 2024-01-01T05:00:00           30.0
   7 │ 2024-01-01T06:00:00           30.0
   8 │ 2024-01-01T07:00:00           30.0
  ⋮  │          ⋮                 ⋮
  18 │ 2024-01-01T17:00:00           30.0
  19 │ 2024-01-01T18:00:00           30.0
  20 │ 2024-01-01T19:00:00           30.0
  21 │ 2024-01-01T20:00:00           30.0
  22 │ 2024-01-01T21:00:00           30.0
  23 │ 2024-01-01T22:00:00           30.0
  24 │ 2024-01-01T23:00:00           30.0
                            9 rows omitted

and the head level of the reservoir:

julia> hydraulic_head = read_variable(
           res,
           "HydroReservoirHeadVariable__HydroReservoir";
           table_format = TableFormat.WIDE,
       )24×2 DataFrame
 Row  DateTime             Water_Reservoir 
      DateTime             Float64?        
─────┼──────────────────────────────────────
   1 │ 2024-01-01T00:00:00          499.718
   2 │ 2024-01-01T01:00:00          499.486
   3 │ 2024-01-01T02:00:00          499.255
   4 │ 2024-01-01T03:00:00          499.023
   5 │ 2024-01-01T04:00:00          498.791
   6 │ 2024-01-01T05:00:00          498.559
   7 │ 2024-01-01T06:00:00          498.327
   8 │ 2024-01-01T07:00:00          498.096
  ⋮  │          ⋮                  ⋮
  18 │ 2024-01-01T17:00:00          495.778
  19 │ 2024-01-01T18:00:00          495.546
  20 │ 2024-01-01T19:00:00          495.314
  21 │ 2024-01-01T20:00:00          495.082
  22 │ 2024-01-01T21:00:00          494.85
  23 │ 2024-01-01T22:00:00          494.619
  24 │ 2024-01-01T23:00:00          494.387
                              9 rows omitted

Note that since the water outflow limit of the turbine is limited on 30 m³/s, the optimal solution decides to flow as much water as possible producing power around 190 MW with that flow and hydraulic head.