Operation Problem with HydroPowerSimulations.jl using Market Bid Cost

To follow along, you can download this tutorial as a Julia script (.jl) or Jupyter notebook (.ipynb).

Load packages

using PowerSystems
using PowerSimulations
using HydroPowerSimulations
using PowerSystemCaseBuilder
using HiGHS ## solver
using TimeSeries
using Dates

Data

sys = build_system(PSITestSystems, "c_sys5_hy"; add_single_time_series = true)

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:

brighton = get_component(ThermalStandard, sys, "Brighton")
old_thermal_cost = get_operation_cost(brighton)
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
)
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,
)
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:

existing_ts = get_time_series_array(
    SingleTimeSeries,
    first(get_components(PowerLoad, sys)),
    "max_active_power",
)
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.

fuel_cost_values = rand(length(tstamps)) .+ 1.0 ## Random fuel cost between 1.0 and 2.0 $/MMBTU
fuel_cost_tarray = TimeArray(tstamps, fuel_cost_values)
fuel_cost_ts = SingleTimeSeries(; name = "fuel_cost", data = fuel_cost_tarray)
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:

existing_ts = get_time_series_array(
    SingleTimeSeries,
    first(get_components(PowerLoad, sys)),
    "max_active_power",
)
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:

hy = get_component(HydroDispatch, sys, "HydroDispatch")
# 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,
)
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:

psd1 = PiecewiseStepData([0.0, 600.0], [5.0])
psd2 = PiecewiseStepData([0.0, 300.0, 600.0], [10.0, 20.0])
psd3 = PiecewiseStepData([0.0, 600.0], [500.0])

# 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])
mbid_tarray = TimeArray(tstamps, total_step_data)
ts_mbid = SingleTimeSeries(; name = "variable_cost", data = mbid_tarray)

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.

zero_input = zeros(length(tstamps))
zero_tarray = TimeArray(tstamps, zero_input)
ts_zero = SingleTimeSeries(; name = "variable_cost_initial_input", data = zero_tarray)
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

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

And create the necessary templates for the system:

template = ProblemTemplate(
    NetworkModel(
        CopperPlatePowerModel;
        use_slacks = true,
        duals = [CopperPlateBalanceConstraint],
    ),
)
set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment)
set_device_model!(template, HydroDispatch, HydroDispatchRunOfRiver)
set_device_model!(template, PowerLoad, StaticPowerLoad)

And then running the decision model:

model = DecisionModel(
    template,
    sys;
    name = "UC_MBCost",
    optimizer = optimizer_with_attributes(
        HiGHS.Optimizer,
    ),
    store_variable_names = true,
    optimizer_solve_log_print = false,
)
build!(model; output_dir = mktempdir())
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.

res = OptimizationProblemResults(model)
hy_p = read_variable(
    res,
    "ActivePowerVariable__HydroDispatch";
    table_format = TableFormat.WIDE,
);
show(hy_p; allrows = true)

dual_price =
    read_dual(res, "CopperPlateBalanceConstraint__System"; table_format = TableFormat.WIDE);
show(dual_price; 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       401.161
  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.024×2 DataFrame
 Row │ DateTime             4        
     │ DateTime             Float64? 
─────┼───────────────────────────────
   1 │ 2024-01-01T00:00:00   1500.0
   2 │ 2024-01-01T01:00:00   1388.06
   3 │ 2024-01-01T02:00:00   1255.45
   4 │ 2024-01-01T03:00:00   1218.87
   5 │ 2024-01-01T04:00:00    957.31
   6 │ 2024-01-01T05:00:00   1020.07
   7 │ 2024-01-01T06:00:00   1400.0
   8 │ 2024-01-01T07:00:00   1083.33
   9 │ 2024-01-01T08:00:00   1500.0
  10 │ 2024-01-01T09:00:00   1543.14
  11 │ 2024-01-01T10:00:00   1379.43
  12 │ 2024-01-01T11:00:00   1053.59
  13 │ 2024-01-01T12:00:00   1106.65
  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