GitHub

Operation Problem with HydroPowerSimulations.jl using Market Bid Cost

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 HiGHS # solver
julia> using TimeSeries
julia> using Dates

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"; add_single_time_series = true)┌ Info: Building new system c_sys5_hy 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/6446c7fb15cee0aebc8ad32e3f600bcdd6f06aab49ede9c3c454114303f80311/c_sys5_hy_time_series_storage.h5.
[ Info: Serialized System to /home/runner/.julia/packages/PowerSystemCaseBuilder/zW01F/data/serialized_system/6446c7fb15cee0aebc8ad32e3f600bcdd6f06aab49ede9c3c454114303f80311/c_sys5_hy.json
[ Info: Serialized System metadata to /home/runner/.julia/packages/PowerSystemCaseBuilder/zW01F/data/serialized_system/6446c7fb15cee0aebc8ad32e3f600bcdd6f06aab49ede9c3c454114303f80311/c_sys5_hy_metadata.json
              System
┌───────────────────┬─────────────┐
│ Property          │ Value       │
├───────────────────┼─────────────┤
│ Name              │             │
│ Description       │             │
│ System Units Base │ SYSTEM_BASE │
│ Base Power        │ 100.0       │
│ Base Frequency    │ 60.0        │
│ Num Components    │ 26          │
└───────────────────┴─────────────┘

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

                            StaticTimeSeries Summary
┌───────────────┬────────────────┬──────────────────┬──────────────────┬────────
│ owner_type    │ owner_category │ name             │ time_series_type │ initi ⋯
│ String        │ String         │ String           │ String           │ Strin ⋯
├───────────────┼────────────────┼──────────────────┼──────────────────┼────────
│ HydroDispatch │ Component      │ max_active_power │ SingleTimeSeries │ 2024- ⋯
│ PowerLoad     │ Component      │ max_active_power │ SingleTimeSeries │ 2024- ⋯
└───────────────┴────────────────┴──────────────────┴──────────────────┴────────
                                                               4 columns omitted

Add a time-varying fuel cost for a thermal unit

We will modify the cheapest unit Brighton to have time-varying fuel cost. First we add a PowerSystems.FuelCurve:

julia> brighton = get_component(ThermalStandard, sys, "Brighton")ThermalStandard: Brighton:
   name: Brighton
   available: true
   status: true
   bus: ACBus: nodeE
   active_power: 6.0
   reactive_power: 1.5
   rating: 0.75
   active_power_limits: (min = 0.0, max = 6.0)
   reactive_power_limits: (min = -4.5, max = 4.5)
   ramp_limits: (up = 0.11249999999999999, down = 0.11249999999999999)
   operation_cost: PowerSystems.ThermalGenerationCost composed of variable: InfrastructureSystems.CostCurve{InfrastructureSystems.LinearCurve}
   base_power: 100.0
   time_limits: (up = 5.0, down = 3.0)
   must_run: false
   prime_mover_type: PowerSystems.PrimeMoversModule.PrimeMovers.ST = 20
   fuel: PowerSystems.ThermalFuelsModule.ThermalFuels.COAL = 1
   services: 0-element Vector{PowerSystems.Service}
   time_at_status: 10000.0
   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> old_thermal_cost = get_operation_cost(brighton)PowerSystems.ThermalGenerationCost:
  variable: CostCurve:
    value_curve: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 10.0 x + 0.0
    power_units: InfrastructureSystems.UnitSystemModule.UnitSystem.NATURAL_UNITS = 2
    vom_cost: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0
  fixed: 0.0
  start_up: 1.5
  shut_down: 0.75
julia> new_fuel_curve = FuelCurve(;
           value_curve = LinearCurve(8.0), # Typical plant of 8 MMBTU/MWh heat rate. Piecewise heat rates can be used if needed.
           power_units = UnitSystem.NATURAL_UNITS,
           fuel_cost = 1.0, # $/MMBTU default fuel cost to start
       )FuelCurve:
  value_curve: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 8.0 x + 0.0
  power_units: InfrastructureSystems.UnitSystemModule.UnitSystem.NATURAL_UNITS = 2
  fuel_cost: 1.0
  startup_fuel_offtake: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0
  vom_cost: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0
julia> new_thermal_cost = ThermalGenerationCost(;
           variable = new_fuel_curve,
           fixed = old_thermal_cost.fixed,
           start_up = old_thermal_cost.start_up,
           shut_down = old_thermal_cost.shut_down,
       )PowerSystems.ThermalGenerationCost:
  variable: FuelCurve:
    value_curve: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 8.0 x + 0.0
    power_units: InfrastructureSystems.UnitSystemModule.UnitSystem.NATURAL_UNITS = 2
    fuel_cost: 1.0
    startup_fuel_offtake: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0
    vom_cost: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0
  fixed: 0.0
  start_up: 1.5
  shut_down: 0.75
julia> set_operation_cost!(brighton, new_thermal_cost)PowerSystems.ThermalGenerationCost:
  variable: FuelCurve:
    value_curve: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 8.0 x + 0.0
    power_units: InfrastructureSystems.UnitSystemModule.UnitSystem.NATURAL_UNITS = 2
    fuel_cost: 1.0
    startup_fuel_offtake: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0
    vom_cost: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0
  fixed: 0.0
  start_up: 1.5
  shut_down: 0.75

Now we create a timeseries with random fuel prices. We first grab an existing timeseries:

julia> existing_ts = get_time_series_array(
           SingleTimeSeries,
           first(get_components(PowerLoad, sys)),
           "max_active_power",
       )48×1 TimeSeries.TimeArray{Float64, 1, Dates.DateTime, Vector{Float64}} 2024-01-01T00:00:00 to 2024-01-02T23:00:00
┌─────────────────────┬─────────┐
│                     │ A       │
├─────────────────────┼─────────┤
│ 2024-01-01T00:00:00 │ 2.37819 │
│ 2024-01-01T01:00:00 │  2.1696 │
│ 2024-01-01T02:00:00 │ 2.13286 │
│ 2024-01-01T03:00:00 │ 2.03302 │
│ 2024-01-01T04:00:00 │ 2.00475 │
│ 2024-01-01T05:00:00 │ 2.01501 │
│ 2024-01-01T06:00:00 │ 2.06283 │
│ 2024-01-01T07:00:00 │ 2.13546 │
│          ⋮          │    ⋮    │
│ 2024-01-02T17:00:00 │ 2.93332 │
│ 2024-01-02T18:00:00 │ 3.03353 │
│ 2024-01-02T19:00:00 │ 3.10266 │
│ 2024-01-02T20:00:00 │ 3.20308 │
│ 2024-01-02T21:00:00 │ 3.05404 │
│ 2024-01-02T22:00:00 │ 2.76541 │
│ 2024-01-02T23:00:00 │ 2.65748 │
└─────────────────────┴─────────┘
                  33 rows omitted
julia> tstamps = timestamp(existing_ts)48-element Vector{Dates.DateTime}:
 2024-01-01T00:00:00
 2024-01-01T01:00:00
 2024-01-01T02:00:00
 2024-01-01T03:00:00
 2024-01-01T04:00:00
 2024-01-01T05:00:00
 2024-01-01T06:00:00
 2024-01-01T07:00:00
 2024-01-01T08:00:00
 2024-01-01T09:00:00
 ⋮
 2024-01-02T15:00:00
 2024-01-02T16:00:00
 2024-01-02T17:00:00
 2024-01-02T18:00:00
 2024-01-02T19:00:00
 2024-01-02T20:00:00
 2024-01-02T21:00:00
 2024-01-02T22:00:00
 2024-01-02T23:00:00

And add the timeseries with the set_fuel_cost! method.

julia> fuel_cost_values = rand(length(tstamps)) .+ 1.0 # Random fuel cost between 1.0 and 2.0 $/MMBTU48-element Vector{Float64}:
 1.3785953776354374
 1.2107644960753656
 1.9337920106658992
 1.336390774754887
 1.744292027840621
 1.0451438385329488
 1.2326450914952156
 1.4569684803092529
 1.8636071171491397
 1.5341504408307534
 ⋮
 1.8209072626779816
 1.7806516529248064
 1.6198708324713829
 1.2525130887454519
 1.2497307026440283
 1.4667076142303146
 1.4422511094309123
 1.2878671209062524
 1.4965707685292138
julia> fuel_cost_tarray = TimeArray(tstamps, fuel_cost_values)48×1 TimeSeries.TimeArray{Float64, 1, Dates.DateTime, Vector{Float64}} 2024-01-01T00:00:00 to 2024-01-02T23:00:00
┌─────────────────────┬─────────┐
│                     │ A       │
├─────────────────────┼─────────┤
│ 2024-01-01T00:00:00 │  1.3786 │
│ 2024-01-01T01:00:00 │ 1.21076 │
│ 2024-01-01T02:00:00 │ 1.93379 │
│ 2024-01-01T03:00:00 │ 1.33639 │
│ 2024-01-01T04:00:00 │ 1.74429 │
│ 2024-01-01T05:00:00 │ 1.04514 │
│ 2024-01-01T06:00:00 │ 1.23265 │
│ 2024-01-01T07:00:00 │ 1.45697 │
│          ⋮          │    ⋮    │
│ 2024-01-02T17:00:00 │ 1.61987 │
│ 2024-01-02T18:00:00 │ 1.25251 │
│ 2024-01-02T19:00:00 │ 1.24973 │
│ 2024-01-02T20:00:00 │ 1.46671 │
│ 2024-01-02T21:00:00 │ 1.44225 │
│ 2024-01-02T22:00:00 │ 1.28787 │
│ 2024-01-02T23:00:00 │ 1.49657 │
└─────────────────────┴─────────┘
                  33 rows omitted
julia> fuel_cost_ts = SingleTimeSeries(; name = "fuel_cost", data = fuel_cost_tarray)InfrastructureSystems.SingleTimeSeries("fuel_cost", 48×1 TimeSeries.TimeArray{Float64, 1, Dates.DateTime, Vector{Float64}} 2024-01-01T00:00:00 to 2024-01-02T23:00:00, Dates.Millisecond(3600000), nothing, InfrastructureSystems.InfrastructureSystemsInternal(Base.UUID("55fc02f5-0027-4851-9a36-ab642f3a661f"), nothing, nothing, nothing))
julia> set_fuel_cost!(sys, brighton, fuel_cost_ts)FuelCurve:
  value_curve: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 8.0 x + 0.0
  power_units: InfrastructureSystems.UnitSystemModule.UnitSystem.NATURAL_UNITS = 2
  fuel_cost: InfrastructureSystems.StaticTimeSeriesKey(InfrastructureSystems.SingleTimeSeries, "fuel_cost", Dates.DateTime("2024-01-01T00:00:00"), Dates.Millisecond(3600000), 48, Dict{String, Any}())
  startup_fuel_offtake: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0
  vom_cost: InfrastructureSystems.LinearCurve (a type of InfrastructureSystems.InputOutputCurve) where function is: f(x) = 0.0 x + 0.0

Add a market bid cost to the hydro unit

We again grab the timestamps from an existing time series:

julia> existing_ts = get_time_series_array(
           SingleTimeSeries,
           first(get_components(PowerLoad, sys)),
           "max_active_power",
       )48×1 TimeSeries.TimeArray{Float64, 1, Dates.DateTime, Vector{Float64}} 2024-01-01T00:00:00 to 2024-01-02T23:00:00
┌─────────────────────┬─────────┐
│                     │ A       │
├─────────────────────┼─────────┤
│ 2024-01-01T00:00:00 │ 2.37819 │
│ 2024-01-01T01:00:00 │  2.1696 │
│ 2024-01-01T02:00:00 │ 2.13286 │
│ 2024-01-01T03:00:00 │ 2.03302 │
│ 2024-01-01T04:00:00 │ 2.00475 │
│ 2024-01-01T05:00:00 │ 2.01501 │
│ 2024-01-01T06:00:00 │ 2.06283 │
│ 2024-01-01T07:00:00 │ 2.13546 │
│          ⋮          │    ⋮    │
│ 2024-01-02T17:00:00 │ 2.93332 │
│ 2024-01-02T18:00:00 │ 3.03353 │
│ 2024-01-02T19:00:00 │ 3.10266 │
│ 2024-01-02T20:00:00 │ 3.20308 │
│ 2024-01-02T21:00:00 │ 3.05404 │
│ 2024-01-02T22:00:00 │ 2.76541 │
│ 2024-01-02T23:00:00 │ 2.65748 │
└─────────────────────┴─────────┘
                  33 rows omitted
julia> tstamps = timestamp(existing_ts)48-element Vector{Dates.DateTime}:
 2024-01-01T00:00:00
 2024-01-01T01:00:00
 2024-01-01T02:00:00
 2024-01-01T03:00:00
 2024-01-01T04:00:00
 2024-01-01T05:00:00
 2024-01-01T06:00:00
 2024-01-01T07:00:00
 2024-01-01T08:00:00
 2024-01-01T09:00:00
 ⋮
 2024-01-02T15:00:00
 2024-01-02T16:00:00
 2024-01-02T17:00:00
 2024-01-02T18:00:00
 2024-01-02T19:00:00
 2024-01-02T20:00:00
 2024-01-02T21:00:00
 2024-01-02T22:00:00
 2024-01-02T23:00:00

And we add an empty market bid cost to the hydro:

julia> hy = get_component(HydroDispatch, sys, "HydroDispatch")HydroDispatch: HydroDispatch:
   name: HydroDispatch
   available: true
   bus: ACBus: nodeB
   active_power: 0.0
   reactive_power: 0.0
   rating: 6.0
   prime_mover_type: PowerSystems.PrimeMoversModule.PrimeMovers.HY = 16
   active_power_limits: (min = 0.0, max = 6.0)
   reactive_power_limits: (min = 0.0, max = 6.0)
   ramp_limits: nothing
   time_limits: nothing
   base_power: 100.0
   status: false
   time_at_status: 10000.0
   operation_cost: PowerSystems.HydroGenerationCost composed of variable: InfrastructureSystems.CostCurve{InfrastructureSystems.LinearCurve}
   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: true
julia> # Create an empty market bid and set it
       hydro_cost = MarketBidCost(;
           no_load_cost = 0.0,
           start_up = (hot = 0.0, warm = 0.0, cold = 0.0),
           shut_down = 0.0,
       )PowerSystems.MarketBidCost:
  no_load_cost: 0.0
  start_up: (hot = 0.0, warm = 0.0, cold = 0.0)
  shut_down: 0.0
  incremental_offer_curves: nothing
  decremental_offer_curves: nothing
  incremental_initial_input: nothing
  decremental_initial_input: nothing
  ancillary_service_offers: PowerSystems.Service[]
julia> set_operation_cost!(hy, hydro_cost)PowerSystems.MarketBidCost:
  no_load_cost: 0.0
  start_up: (hot = 0.0, warm = 0.0, cold = 0.0)
  shut_down: 0.0
  incremental_offer_curves: nothing
  decremental_offer_curves: nothing
  incremental_initial_input: nothing
  decremental_initial_input: nothing
  ancillary_service_offers: PowerSystems.Service[]

Now we create a time-varying piecewise linear bid cost:

julia> psd1 = PiecewiseStepData([0.0, 600.0], [5.0])InfrastructureSystems.PiecewiseStepData representing step (piecewise constant) function f(x) =
  5.0 for x in [0.0, 600.0)
julia> psd2 = PiecewiseStepData([0.0, 300.0, 600.0], [10.0, 20.0])InfrastructureSystems.PiecewiseStepData representing step (piecewise constant) function f(x) =
  10.0 for x in [0.0, 300.0)
  20.0 for x in [300.0, 600.0)
julia> psd3 = PiecewiseStepData([0.0, 600.0], [500.0])InfrastructureSystems.PiecewiseStepData representing step (piecewise constant) function f(x) =
  500.0 for x in [0.0, 600.0)
julia> # Cheap the first 10 hours, moderate next 4 hours, expensive last 34 hours
       total_step_data = vcat([psd1 for x in 1:10], [psd2 for x in 1:4], [psd3 for x in 1:34])48-element Vector{InfrastructureSystems.PiecewiseStepData}:
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [5.0])
 ⋮
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
 InfrastructureSystems.PiecewiseStepData([0.0, 600.0], [500.0])
julia> mbid_tarray = TimeArray(tstamps, total_step_data)48×1 TimeSeries.TimeArray{InfrastructureSystems.PiecewiseStepData, 1, Dates.DateTime, Vector{InfrastructureSystems.PiecewiseStepData}} 2024-01-01T00:00:00 to 2024-01-02T23:00:00
┌─────────────────────┬──────────────────────────────────────────┐
│                     │ A                                        │
├─────────────────────┼──────────────────────────────────────────┤
│ 2024-01-01T00:00:00 │   PiecewiseStepData([0.0, 600.0], [5.0]) │
│ 2024-01-01T01:00:00 │   PiecewiseStepData([0.0, 600.0], [5.0]) │
│ 2024-01-01T02:00:00 │   PiecewiseStepData([0.0, 600.0], [5.0]) │
│ 2024-01-01T03:00:00 │   PiecewiseStepData([0.0, 600.0], [5.0]) │
│ 2024-01-01T04:00:00 │   PiecewiseStepData([0.0, 600.0], [5.0]) │
│ 2024-01-01T05:00:00 │   PiecewiseStepData([0.0, 600.0], [5.0]) │
│ 2024-01-01T06:00:00 │   PiecewiseStepData([0.0, 600.0], [5.0]) │
│ 2024-01-01T07:00:00 │   PiecewiseStepData([0.0, 600.0], [5.0]) │
│          ⋮          │                    ⋮                     │
│ 2024-01-02T17:00:00 │ PiecewiseStepData([0.0, 600.0], [500.0]) │
│ 2024-01-02T18:00:00 │ PiecewiseStepData([0.0, 600.0], [500.0]) │
│ 2024-01-02T19:00:00 │ PiecewiseStepData([0.0, 600.0], [500.0]) │
│ 2024-01-02T20:00:00 │ PiecewiseStepData([0.0, 600.0], [500.0]) │
│ 2024-01-02T21:00:00 │ PiecewiseStepData([0.0, 600.0], [500.0]) │
│ 2024-01-02T22:00:00 │ PiecewiseStepData([0.0, 600.0], [500.0]) │
│ 2024-01-02T23:00:00 │ PiecewiseStepData([0.0, 600.0], [500.0]) │
└─────────────────────┴──────────────────────────────────────────┘
                                                   33 rows omitted
julia> ts_mbid = SingleTimeSeries(; name = "variable_cost", data = mbid_tarray)InfrastructureSystems.SingleTimeSeries("variable_cost", 48×1 TimeSeries.TimeArray{InfrastructureSystems.PiecewiseStepData, 1, Dates.DateTime, Vector{InfrastructureSystems.PiecewiseStepData}} 2024-01-01T00:00:00 to 2024-01-02T23:00:00, Dates.Millisecond(3600000), nothing, InfrastructureSystems.InfrastructureSystemsInternal(Base.UUID("d1fc9c94-1e67-4beb-85c9-5c2b29e19943"), nothing, nothing, nothing))
julia> set_variable_cost!(sys, hy, ts_mbid, UnitSystem.NATURAL_UNITS)

It is also needed to create the initial input time series for market bid. That is the cost at 0 power at each time step. We will use zero for this example.

julia> zero_input = zeros(length(tstamps))48-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
julia> zero_tarray = TimeArray(tstamps, zero_input)48×1 TimeSeries.TimeArray{Float64, 1, Dates.DateTime, Vector{Float64}} 2024-01-01T00:00:00 to 2024-01-02T23:00:00
┌─────────────────────┬─────┐
│                     │ A   │
├─────────────────────┼─────┤
│ 2024-01-01T00:00:00 │ 0.0 │
│ 2024-01-01T01:00:00 │ 0.0 │
│ 2024-01-01T02:00:00 │ 0.0 │
│ 2024-01-01T03:00:00 │ 0.0 │
│ 2024-01-01T04:00:00 │ 0.0 │
│ 2024-01-01T05:00:00 │ 0.0 │
│ 2024-01-01T06:00:00 │ 0.0 │
│ 2024-01-01T07:00:00 │ 0.0 │
│          ⋮          │  ⋮  │
│ 2024-01-02T17:00:00 │ 0.0 │
│ 2024-01-02T18:00:00 │ 0.0 │
│ 2024-01-02T19:00:00 │ 0.0 │
│ 2024-01-02T20:00:00 │ 0.0 │
│ 2024-01-02T21:00:00 │ 0.0 │
│ 2024-01-02T22:00:00 │ 0.0 │
│ 2024-01-02T23:00:00 │ 0.0 │
└─────────────────────┴─────┘
              33 rows omitted
julia> ts_zero = SingleTimeSeries(; name = "variable_cost_initial_input", data = zero_tarray)InfrastructureSystems.SingleTimeSeries("variable_cost_initial_input", 48×1 TimeSeries.TimeArray{Float64, 1, Dates.DateTime, Vector{Float64}} 2024-01-01T00:00:00 to 2024-01-02T23:00:00, Dates.Millisecond(3600000), nothing, InfrastructureSystems.InfrastructureSystemsInternal(Base.UUID("ef68b848-20b2-461d-9d84-a3af268118a2"), nothing, nothing, nothing))
julia> set_incremental_initial_input!(sys, hy, ts_zero)InfrastructureSystems.StaticTimeSeriesKey(InfrastructureSystems.SingleTimeSeries, "variable_cost_initial_input", Dates.DateTime("2024-01-01T00:00:00"), Dates.Millisecond(3600000), 48, Dict{String, Any}())

Running the single-stage problem

We first transform the single time series to a 24 hour forecast

julia> transform_single_time_series!(sys, Hour(24), Hour(24))

And create the necessary templates for the system:

julia> template = ProblemTemplate(
           NetworkModel(
               CopperPlatePowerModel;
               use_slacks = true,
               duals = [CopperPlateBalanceConstraint],
           ),
       )                            Network Model
┌────────────────────┬───────────────────────────────────────────────┐
│ Network Model      │ PowerSimulations.CopperPlatePowerModel        │
│ Slacks             │ true                                          │
│ PTDF               │ false                                         │
│ Duals              │ PowerSimulations.CopperPlateBalanceConstraint │
│ HVDC Network Model │ None                                          │
└────────────────────┴───────────────────────────────────────────────┘

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

And then running the decision model:

julia> model = DecisionModel(
           template,
           sys;
           name = "UC_MBCost",
           optimizer = optimizer_with_attributes(
               HiGHS.Optimizer,
           ),
           store_variable_names = true,
           optimizer_solve_log_print = false,
       )[ Info: Overriding time_series_cache_size because time series is stored in memory
                            Network Model
┌────────────────────┬───────────────────────────────────────────────┐
│ Network Model      │ PowerSimulations.CopperPlatePowerModel        │
│ Slacks             │ true                                          │
│ PTDF               │ false                                         │
│ Duals              │ PowerSimulations.CopperPlateBalanceConstraint │
│ HVDC Network Model │ None                                          │
└────────────────────┴───────────────────────────────────────────────┘

                                 Device Models
┌──────────────────────────────┬─────────────────────────────────────────────┬──
│ Device Type                  │ Formulation                                 │ ⋯
├──────────────────────────────┼─────────────────────────────────────────────┼──
│ PowerSystems.ThermalStandard │ PowerSimulations.ThermalBasicUnitCommitment │ ⋯
│ PowerSystems.PowerLoad       │ PowerSimulations.StaticPowerLoad            │ ⋯
│ PowerSystems.HydroDispatch   │ HydroDispatchRunOfRiver                     │ ⋯
└──────────────────────────────┴─────────────────────────────────────────────┴──
                                                                1 column omitted
julia> build!(model; output_dir = mktempdir())InfrastructureSystems.Optimization.ModelBuildStatusModule.ModelBuildStatus.BUILT = 0
julia> solve!(model)InfrastructureSystems.Simulation.RunStatusModule.RunStatus.SUCCESSFULLY_FINALIZED = 0

And exploring results we confirm that the hydro is not dispatched when is more expensive, while the dual of the CopperPlate constraint showcase that the system become more expensive at those hours.

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

PowerSimulations Problem Auxiliary variables Results
┌──────────────────────────────────┐
│ TimeDurationOff__ThermalStandard │
│ TimeDurationOn__ThermalStandard  │
│ HydroEnergyOutput__HydroDispatch │
└──────────────────────────────────┘

 PowerSimulations Problem Expressions Results
┌────────────────────────────────────────────┐
│ ProductionCostExpression__ThermalStandard  │
│ ActivePowerBalance__System                 │
│ ProductionCostExpression__HydroDispatch    │
│ FuelConsumptionExpression__ThermalStandard │
└────────────────────────────────────────────┘

 PowerSimulations Problem Duals Results
┌──────────────────────────────────────┐
│ CopperPlateBalanceConstraint__System │
└──────────────────────────────────────┘

          PowerSimulations Problem Parameters Results
┌──────────────────────────────────────────────────────────────┐
│ ActivePowerTimeSeriesParameter__HydroDispatch                │
│ IncrementalPiecewiseLinearSlopeParameter__HydroDispatch      │
│ FuelCostParameter__ThermalStandard                           │
│ ActivePowerTimeSeriesParameter__PowerLoad                    │
│ IncrementalPiecewiseLinearBreakpointParameter__HydroDispatch │
└──────────────────────────────────────────────────────────────┘

PowerSimulations Problem Variables Results
┌──────────────────────────────────────┐
│ SystemBalanceSlackUp__System         │
│ ActivePowerVariable__ThermalStandard │
│ StartVariable__ThermalStandard       │
│ OnVariable__ThermalStandard          │
│ SystemBalanceSlackDown__System       │
│ StopVariable__ThermalStandard        │
│ ActivePowerVariable__HydroDispatch   │
└──────────────────────────────────────┘
julia> hy_p = read_variable(
           res,
           "ActivePowerVariable__HydroDispatch";
           table_format = TableFormat.WIDE,
       );
julia> show(hy_p; allrows = true)24×2 DataFrame
 Row │ DateTime             HydroDispatch 
     │ DateTime             Float64?      
─────┼────────────────────────────────────
   1 │ 2024-01-01T00:00:00       188.58
   2 │ 2024-01-01T01:00:00       232.01
   3 │ 2024-01-01T02:00:00       137.149
   4 │ 2024-01-01T03:00:00       136.006
   5 │ 2024-01-01T04:00:00       133.72
   6 │ 2024-01-01T05:00:00        77.718
   7 │ 2024-01-01T06:00:00        86.8608
   8 │ 2024-01-01T07:00:00       219.439
   9 │ 2024-01-01T08:00:00       124.577
  10 │ 2024-01-01T09:00:00       373.731
  11 │ 2024-01-01T10:00:00       402.304
  12 │ 2024-01-01T11:00:00       405.733
  13 │ 2024-01-01T12:00:00       285.644
  14 │ 2024-01-01T13:00:00       244.583
  15 │ 2024-01-01T14:00:00         0.0
  16 │ 2024-01-01T15:00:00         0.0
  17 │ 2024-01-01T16:00:00         0.0
  18 │ 2024-01-01T17:00:00         0.0
  19 │ 2024-01-01T18:00:00         0.0
  20 │ 2024-01-01T19:00:00         0.0
  21 │ 2024-01-01T20:00:00         0.0
  22 │ 2024-01-01T21:00:00         0.0
  23 │ 2024-01-01T22:00:00         0.0
  24 │ 2024-01-01T23:00:00         0.0
julia> dual_price =
           read_dual(res, "CopperPlateBalanceConstraint__System"; table_format = TableFormat.WIDE);
julia> show(dual_price; allrows = true)24×2 DataFrame
 Row │ DateTime             4        
     │ DateTime             Float64? 
─────┼───────────────────────────────
   1 │ 2024-01-01T00:00:00  1500.0
   2 │ 2024-01-01T01:00:00   968.612
   3 │ 2024-01-01T02:00:00  1547.03
   4 │ 2024-01-01T03:00:00  1069.11
   5 │ 2024-01-01T04:00:00  1395.43
   6 │ 2024-01-01T05:00:00   836.115
   7 │ 2024-01-01T06:00:00  1400.0
   8 │ 2024-01-01T07:00:00  1165.57
   9 │ 2024-01-01T08:00:00  1500.0
  10 │ 2024-01-01T09:00:00  1227.32
  11 │ 2024-01-01T10:00:00  1013.29
  12 │ 2024-01-01T11:00:00  1376.34
  13 │ 2024-01-01T12:00:00  1000.0
  14 │ 2024-01-01T13:00:00  1400.0
  15 │ 2024-01-01T14:00:00  3000.0
  16 │ 2024-01-01T15:00:00  3000.0
  17 │ 2024-01-01T16:00:00  3000.0
  18 │ 2024-01-01T17:00:00  3000.0
  19 │ 2024-01-01T18:00:00  3000.0
  20 │ 2024-01-01T19:00:00  3000.0
  21 │ 2024-01-01T20:00:00  3000.0
  22 │ 2024-01-01T21:00:00  3000.0
  23 │ 2024-01-01T22:00:00  3000.0
  24 │ 2024-01-01T23:00:00  1500.0