Sequential Simulations with PowerSimulations.jl

Originally Contributed by: Clayton Barrows

Introduction

PowerSimulations.jl supports simulations that consist of sequential optimization problems where results from previous problems inform subsequent problems in a variety of ways. This example demonstrates some of these capabilities to represent electricity market clearing. This example is intended to be an extension of the OperationsProblem tutorial.

Load Packages

using PowerSystems
using PowerSimulations
using HydroPowerSimulations
const PSI = PowerSimulations
using PowerSystemCaseBuilder
using Dates
using HiGHS #solver

Optimizer

It's most convenient to define an optimizer instance upfront and pass it into the DecisionModel constructor. For this example, we can use the free HiGHS solver with a relatively relaxed MIP gap (ratioGap) setting to improve speed.

solver = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.5)
MathOptInterface.OptimizerWithAttributes(HiGHS.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("mip_rel_gap") => 0.5])

Hourly day-ahead system

First, we'll create a System with hourly data to represent day-ahead forecasted wind, solar, and load profiles:

sys_DA = build_system(PSISystems, "modified_RTS_GMLC_DA_sys"; skip_serialization = true)
System
Property Value
Name
Description
System Units Base SYSTEM_BASE
Base Power 100.0
Base Frequency 60.0
Num Components 501
Static Components
Type Count
ACBus 73
Arc 109
Area 3
FixedAdmittance 3
HydroDispatch 1
Line 105
LoadZone 21
PowerLoad 51
RenewableDispatch 29
RenewableNonDispatch 31
TapTransformer 15
ThermalStandard 54
TwoTerminalHVDCLine 1
VariableReserve{ReserveDown} 1
VariableReserve{ReserveUp} 4
StaticTimeSeries Summary
owner_type owner_category time_series_type initial_timestamp resolution count time_step_count
String String String String Dates.CompoundPeriod Int64 Int64
Area Component SingleTimeSeries 2020-01-01T00:00:00 1 hour 3 8784
FixedAdmittance Component SingleTimeSeries 2020-01-01T00:00:00 1 hour 3 8784
HydroDispatch Component SingleTimeSeries 2020-01-01T00:00:00 1 hour 2 8784
PowerLoad Component SingleTimeSeries 2020-01-01T00:00:00 1 hour 51 8784
RenewableDispatch Component SingleTimeSeries 2020-01-01T00:00:00 1 hour 29 8784
RenewableNonDispatch Component SingleTimeSeries 2020-01-01T00:00:00 1 hour 31 8784
VariableReserve Component SingleTimeSeries 2020-01-01T00:00:00 1 hour 5 8784
Forecast Summary
owner_type owner_category time_series_type initial_timestamp resolution count horizon interval window_count
String String String String Dates.CompoundPeriod Int64 Dates.CompoundPeriod Dates.CompoundPeriod Int64
Area Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 1 hour 3 2 days 1 day 365
FixedAdmittance Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 1 hour 3 2 days 1 day 365
HydroDispatch Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 1 hour 2 2 days 1 day 365
PowerLoad Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 1 hour 51 2 days 1 day 365
RenewableDispatch Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 1 hour 29 2 days 1 day 365
RenewableNonDispatch Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 1 hour 31 2 days 1 day 365
VariableReserve Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 1 hour 5 2 days 1 day 365

5-Minute system

The RTS data also includes 5-minute resolution time series data. So, we can create another System to represent 15 minute ahead forecasted data for a "real-time" market:

sys_RT = build_system(PSISystems, "modified_RTS_GMLC_RT_sys"; skip_serialization = true)
System
Property Value
Name
Description
System Units Base SYSTEM_BASE
Base Power 100.0
Base Frequency 60.0
Num Components 499
Static Components
Type Count
ACBus 73
Arc 109
Area 1
FixedAdmittance 3
HydroDispatch 1
Line 105
LoadZone 21
PowerLoad 51
RenewableDispatch 29
RenewableNonDispatch 31
TapTransformer 15
ThermalStandard 54
TwoTerminalHVDCLine 1
VariableReserve{ReserveDown} 1
VariableReserve{ReserveUp} 4
StaticTimeSeries Summary
owner_type owner_category time_series_type initial_timestamp resolution count time_step_count
String String String String Dates.CompoundPeriod Int64 Int64
Area Component SingleTimeSeries 2020-01-01T00:00:00 5 minutes 1 105408
FixedAdmittance Component SingleTimeSeries 2020-01-01T00:00:00 5 minutes 3 105408
HydroDispatch Component SingleTimeSeries 2020-01-01T00:00:00 5 minutes 2 105408
PowerLoad Component SingleTimeSeries 2020-01-01T00:00:00 5 minutes 51 105408
RenewableDispatch Component SingleTimeSeries 2020-01-01T00:00:00 5 minutes 29 105408
RenewableNonDispatch Component SingleTimeSeries 2020-01-01T00:00:00 5 minutes 31 105408
VariableReserve Component SingleTimeSeries 2020-01-01T00:00:00 5 minutes 5 105408
Forecast Summary
owner_type owner_category time_series_type initial_timestamp resolution count horizon interval window_count
String String String String Dates.CompoundPeriod Int64 Dates.CompoundPeriod Dates.CompoundPeriod Int64
Area Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 5 minutes 1 1 hour 15 minutes 35133
FixedAdmittance Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 5 minutes 3 1 hour 15 minutes 35133
HydroDispatch Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 5 minutes 2 1 hour 15 minutes 35133
PowerLoad Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 5 minutes 51 1 hour 15 minutes 35133
RenewableDispatch Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 5 minutes 29 1 hour 15 minutes 35133
RenewableNonDispatch Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 5 minutes 31 1 hour 15 minutes 35133
VariableReserve Component DeterministicSingleTimeSeries 2020-01-01T00:00:00 5 minutes 5 1 hour 15 minutes 35133

ProblemTemplates define stages

Sequential simulations in PowerSimulations are created by defining OperationsProblems that represent stages, and how information flows between executions of a stage and between different stages.

Let's start by defining a two stage simulation that might look like a typical day-Ahead and real-time electricity market clearing process.

Day-ahead unit commitment stage

First, we can define a unit commitment template for the day ahead problem. We can use the included UC template, but in this example, we'll replace the ThermalBasicUnitCommitment with the slightly more complex ThermalStandardUnitCommitment for the thermal generators.

template_uc = template_unit_commitment()
set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment)
set_device_model!(template_uc, HydroDispatch, HydroDispatchRunOfRiver)
┌ Warning: Overwriting ThermalStandard existing model
@ PowerSimulations ~/work/PowerSimulations.jl/PowerSimulations.jl/src/core/device_model.jl:116

Define the reference model for the real-time economic dispatch

In addition to the manual specification process demonstrated in the OperationsProblem example, PSI also provides pre-specified templates for some standard problems:

template_ed = template_economic_dispatch(;
    network = NetworkModel(PTDFPowerModel; use_slacks = true),
)
Network Model
Network Model PTDFPowerModel
Slacks true
PTDF false
Duals None
Device Models
Device Type Formulation Slacks
RenewableNonDispatch FixedOutput false
ThermalStandard ThermalBasicDispatch false
PowerLoad StaticPowerLoad false
InterruptiblePowerLoad PowerLoadInterruption false
RenewableDispatch RenewableFullDispatch false
Branch Models
Branch Type Formulation Slacks
Line StaticBranch false
TapTransformer StaticBranch false
Transformer2W StaticBranch false
TwoTerminalHVDCLine HVDCTwoTerminalDispatch false
Service Models
Service Type Formulation Slacks Aggregated Model
VariableReserve{ReserveUp} RangeReserve false true
VariableReserve{ReserveDown} RangeReserve false true

Define the SimulationModels

DecisionModels define the problems that are executed in the simulation. The actual problem will change as the stage gets updated to represent different time periods, but the formulations applied to the components is constant within a stage. In this case, we want to define two stages with the ProblemTemplates and the Systems that we've already created.

models = SimulationModels(;
    decision_models = [
        DecisionModel(template_uc, sys_DA; optimizer = solver, name = "UC"),
        DecisionModel(template_ed, sys_RT; optimizer = solver, name = "ED"),
    ],
)
Decision Models
Model Name Model Type Status Output Directory
UC GenericOpProblem EMPTY nothing
ED GenericOpProblem EMPTY nothing
No Emulator Model Specified

SimulationSequence

Similar to an ProblemTemplate, the SimulationSequence provides a template of how to execute a sequential set of operations problems.

Let's review some of the SimulationSequence arguments.

Chronologies

In PowerSimulations, chronologies define where information is flowing. There are two types of chronologies.

  • inter-stage chronologies: Define how information flows between stages. e.g. day-ahead solutions are used to inform economic dispatch problems
  • intra-stage chronologies: Define how information flows between multiple executions of a single stage. e.g. the dispatch setpoints of the first period of an economic dispatch problem are constrained by the ramping limits from setpoints in the final period of the previous problem.

FeedForward

The definition of exactly what information is passed using the defined chronologies is accomplished with FeedForward. Specifically, FeedForward is used to define what to do with information being passed with an inter-stage chronology. Let's define a FeedForward that affects the semi-continuous range constraints of thermal generators in the economic dispatch problems based on the value of the unit-commitment variables.

feedforward = Dict(
    "ED" => [
        SemiContinuousFeedforward(;
            component_type = ThermalStandard,
            source = OnVariable,
            affected_values = [ActivePowerVariable],
        ),
    ],
)
Dict{String, Vector{SemiContinuousFeedforward}} with 1 entry:
  "ED" => [SemiContinuousFeedforward(VariableKey{OnVariable, ThermalStandard}("…

Sequencing

The stage problem length, look-ahead, and other details surrounding the temporal Sequencing of stages are controlled using the structure of the time series data in the Systems. So, to define a typical day-ahead - real-time sequence:

  • Day ahead problems should represent 48 hours, advancing 24 hours after each execution (24-hour look-ahead)
  • Real time problems should represent 1 hour (12 5-minute periods), advancing 15 min after each execution (15 min look-ahead)

We can adjust the time series data to reflect this structure in each System:

  • transform_single_time_series!(sys_DA, 48, Hour(1))
  • transform_single_time_series!(sys_RT, 12, Minute(15))

Now we can put it all together to define a SimulationSequence

DA_RT_sequence = SimulationSequence(;
    models = models,
    ini_cond_chronology = InterProblemChronology(),
    feedforwards = feedforward,
)
Simulation Sequence
Simulation Step Interval 24 hours
Number of Problems 2
Simulation Problems
Model Name Horizon Interval Executions Per Step
UC 2 days 1 day 1
ED 1 hour 15 minutes 96
Feedforwards
Model Name Feed Forward Type
ED SemiContinuousFeedforward

Simulation

Now, we can build and execute a simulation using the SimulationSequence and Stages that we've defined.

sim = Simulation(;
    name = "rts-test",
    steps = 2,
    models = models,
    sequence = DA_RT_sequence,
    simulation_folder = joinpath(".", "rts-store"),
)
Simulation
Simulation Name rts-test
Build Status EMPTY
Run Status NOT_READY
Initial Time Unset Initial Time
Steps 2
Decision Models
Model Name Model Type Status Output Directory
UC GenericOpProblem EMPTY nothing
ED GenericOpProblem EMPTY nothing
No Emulator Model Specified
Simulation Sequence
Simulation Step Interval 24 hours
Number of Problems 2
Simulation Problems
Model Name Horizon Interval Executions Per Step
UC 2 days 1 day 1
ED 1 hour 15 minutes 96
Feedforwards
Model Name Feed Forward Type
ED SemiContinuousFeedforward

Build simulation

build!(sim)
InfrastructureSystems.Simulation.SimulationBuildStatusModule.SimulationBuildStatus.BUILT = 0

Execute simulation

the following command returns the status of the simulation (0: is proper execution) and stores the results in a set of HDF5 files on disk.

execute!(sim; enable_progress_bar = false)
InfrastructureSystems.Simulation.RunStatusModule.RunStatus.SUCCESSFULLY_FINALIZED = 0

Results

To access the results, we need to load the simulation result metadata and then make requests to the specific data of interest. This allows you to efficiently access the results of interest without overloading resources.

results = SimulationResults(sim);
uc_results = get_decision_problem_results(results, "UC"); # UC stage result metadata
ed_results = get_decision_problem_results(results, "ED"); # ED stage result metadata

Start: 2020-01-01T00:00:00

End: 2020-01-02T23:45:00

Resolution: 5 minutes

ED Problem Expressions Results
ProductionCostExpression__RenewableDispatch
ProductionCostExpression__ThermalStandard
ActivePowerBalance__System
ActivePowerBalance__ACBus
FuelConsumptionExpression__ThermalStandard
ED Problem Parameters Results
OnStatusParameter__ThermalStandard
RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R3
RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R1
ActivePowerTimeSeriesParameter__PowerLoad
ActivePowerTimeSeriesParameter__RenewableNonDispatch
RequirementTimeSeriesParameter__VariableReserve__ReserveDown__Reg_Down
RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R2
ActivePowerTimeSeriesParameter__RenewableDispatch
RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Reg_Up
ED Problem Variables Results
HVDCLosses__TwoTerminalHVDCLine
SystemBalanceSlackUp__System
HVDCFlowDirectionVariable__TwoTerminalHVDCLine
ActivePowerReserveVariable__VariableReserve__ReserveUp__Reg_Up
ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R3
SystemBalanceSlackDown__System
FlowActivePowerFromToVariable__TwoTerminalHVDCLine
ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R1
ActivePowerVariable__ThermalStandard
ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R2
FlowActivePowerToFromVariable__TwoTerminalHVDCLine
FlowActivePowerVariable__Line
ActivePowerVariable__RenewableDispatch
FlowActivePowerVariable__TapTransformer
ActivePowerReserveVariable__VariableReserve__ReserveDown__Reg_Down

We can read all the result variables

read_variables(uc_results)
Dict{String, SortedDict{Any, Any, Base.Order.ForwardOrdering}} with 15 entries:
  "ActivePowerReserveVaria… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×5…
  "FlowActivePowerToFromVa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "StopVariable__ThermalSt… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×5…
  "ActivePowerReserveVaria… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×5…
  "OnVariable__ThermalStan… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×5…
  "ActivePowerVariable__Hy… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "ActivePowerReserveVaria… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×1…
  "StartVariable__ThermalS… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×5…
  "ActivePowerVariable__Th… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×5…
  "FlowActivePowerFromToVa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "ActivePowerReserveVaria… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×1…
  "ActivePowerVariable__Re… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×3…
  "HVDCLosses__TwoTerminal… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "HVDCFlowDirectionVariab… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "ActivePowerReserveVaria… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×1

or all the parameters

read_parameters(uc_results)
Dict{String, SortedDict{Any, Any, Base.Order.ForwardOrdering}} with 9 entries:
  "RequirementTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "ActivePowerTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×5…
  "ActivePowerTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×3…
  "ActivePowerTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×3…
  "RequirementTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "ActivePowerTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "RequirementTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "RequirementTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2…
  "RequirementTimeSeriesPa… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2

We can just list the variable names contained in uc_results:

list_variable_names(uc_results)
15-element Vector{String}:
 "HVDCLosses__TwoTerminalHVDCLine"
 "HVDCFlowDirectionVariable__TwoTerminalHVDCLine"
 "StopVariable__ThermalStandard"
 "ActivePowerReserveVariable__VariableReserve__ReserveUp__Reg_Up"
 "ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R3"
 "FlowActivePowerFromToVariable__TwoTerminalHVDCLine"
 "ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R1"
 "ActivePowerVariable__ThermalStandard"
 "ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R2"
 "FlowActivePowerToFromVariable__TwoTerminalHVDCLine"
 "ActivePowerVariable__HydroDispatch"
 "ActivePowerVariable__RenewableDispatch"
 "OnVariable__ThermalStandard"
 "StartVariable__ThermalStandard"
 "ActivePowerReserveVariable__VariableReserve__ReserveDown__Reg_Down"

and a number of parameters (this pattern also works for aux_variables, expressions, and duals)

list_parameter_names(uc_results)
9-element Vector{String}:
 "RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R3"
 "RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R1"
 "ActivePowerTimeSeriesParameter__PowerLoad"
 "ActivePowerTimeSeriesParameter__RenewableNonDispatch"
 "RequirementTimeSeriesParameter__VariableReserve__ReserveDown__Reg_Down"
 "RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R2"
 "ActivePowerTimeSeriesParameter__HydroDispatch"
 "ActivePowerTimeSeriesParameter__RenewableDispatch"
 "RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Reg_Up"

Now we can read the specific results of interest for a specific problem, time window (optional), and set of variables, duals, or parameters (optional)

Dict([
    v => read_variable(uc_results, v) for v in [
        "ActivePowerVariable__RenewableDispatch",
        "ActivePowerVariable__HydroDispatch",
        "StopVariable__ThermalStandard",
    ]
])
Dict{String, SortedDict{Any, Any, Base.Order.ForwardOrdering}} with 3 entries:
  "StopVariable__ThermalSt… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×5…
  "ActivePowerVariable__Re… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×3…
  "ActivePowerVariable__Hy… => SortedDict(DateTime("2020-01-01T00:00:00")=>48×2

Or if we want the result of just one variable, parameter, or dual (must be defined in the problem definition), we can use:

read_parameter(
    ed_results,
    "ActivePowerTimeSeriesParameter__RenewableNonDispatch";
    initial_time = DateTime("2020-01-01T06:00:00"),
    count = 5,
)
SortedDict{Any, Any, Base.Order.ForwardOrdering} with 5 entries:
  DateTime("2020-01-01T06:00:00") => 12×32 DataFrame…
  DateTime("2020-01-01T06:15:00") => 12×32 DataFrame…
  DateTime("2020-01-01T06:30:00") => 12×32 DataFrame…
  DateTime("2020-01-01T06:45:00") => 12×32 DataFrame…
  DateTime("2020-01-01T07:00:00") => 12×32 DataFrame
Info

note that this returns the results of each execution step in a separate dataframe If you want the realized results (without lookahead periods), you can call read_realized_*:

read_realized_variables(
    uc_results,
    ["ActivePowerVariable__ThermalStandard", "ActivePowerVariable__RenewableDispatch"],
)

Plotting

Take a look at the plotting capabilities in PowerGraphics.jl