Operations problems with PowerSimulations.jl

Originally Contributed by: Clayton Barrows

Introduction

PowerSimulations.jl supports the construction and solution of optimal power system scheduling problems (Operations Problems). Operations problems form the fundamental building blocks for sequential simulations. This example shows how to specify and customize a the mathematics that will be applied to the data with an ProblemTemplate, build and execute an DecisionModel, and access the results.

Load Packages

using PowerSystems
using PowerSimulations
using HydroPowerSimulations
using PowerSystemCaseBuilder
using HiGHS # solver
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 Documentation

sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
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

Define a problem specification with an ProblemTemplate

You can create an empty template with:

template_uc = ProblemTemplate()
Network Model
Network Model CopperPlatePowerModel
Slacks false
PTDF false
Duals None
Device Models
Device Type Formulation Slacks

Now, you can add a DeviceModel for each device type to create an assignment between PowerSystems device types and the subtypes of AbstractDeviceFormulation. PowerSimulations has a variety of different AbstractDeviceFormulation subtypes that can be applied to different PowerSystems device types, each dispatching to different methods for populating optimization problem objectives, variables, and constraints. Documentation on the formulation options for various devices can be found in the formulation library docs

Branch Formulations

Here is an example of relatively standard branch formulations. Other formulations allow for selective enforcement of transmission limits and greater control on transformer settings.

set_device_model!(template_uc, Line, StaticBranch)
set_device_model!(template_uc, Transformer2W, StaticBranch)
set_device_model!(template_uc, TapTransformer, StaticBranch)

Injection Device Formulations

Here we define template entries for all devices that inject or withdraw power on the network. For each device type, we can define a distinct AbstractDeviceFormulation. In this case, we're defining a basic unit commitment model for thermal generators, curtailable renewable generators, and fixed dispatch (net-load reduction) formulations for HydroDispatch and RenewableNonDispatch devices.

set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment)
set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch)
set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
set_device_model!(template_uc, HydroDispatch, HydroDispatchRunOfRiver)
set_device_model!(template_uc, RenewableNonDispatch, FixedOutput)

Service Formulations

We have two VariableReserve types, parameterized by their direction. So, similar to creating DeviceModels, we can create ServiceModels. The primary difference being that DeviceModel objects define how constraints get created, while ServiceModel objects define how constraints get modified.

set_service_model!(template_uc, VariableReserve{ReserveUp}, RangeReserve)
set_service_model!(template_uc, VariableReserve{ReserveDown}, RangeReserve)

Network Formulations

Finally, we can define the transmission network specification that we'd like to model. For simplicity, we'll choose a copper plate formulation. But there are dozens of specifications available through an integration with PowerModels.jl. Note that many formulations will require appropriate data and may be computationally intractable

set_network_model!(template_uc, NetworkModel(CopperPlatePowerModel))

DecisionModel

Now that we have a System and an ProblemTemplate, we can put the two together to create an DecisionModel that we solve.

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])

Build an DecisionModel

The construction of an DecisionModel essentially applies an ProblemTemplate to System data to create a JuMP model.

problem = DecisionModel(template_uc, sys; optimizer = solver, horizon = Hour(24))
build!(problem; output_dir = mktempdir())
InfrastructureSystems.Optimization.ModelBuildStatusModule.ModelBuildStatus.BUILT = 0
Tip

The principal component of the DecisionModel is the JuMP model. But you can serialize to a file using the following command:

serialize_optimization_model(problem, save_path)

Keep in mind that if the setting "storevariablenames" is set to False then the file won't show the model's names.

Solve an DecisionModel

solve!(problem)
InfrastructureSystems.Simulation.RunStatusModule.RunStatus.SUCCESSFULLY_FINALIZED = 0

Results Inspection

PowerSimulations collects the DecisionModel results into a OptimizationProblemResults struct:

res = OptimizationProblemResults(problem)

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

End: 2020-01-01T23:00:00

Resolution: 60 minutes

PowerSimulations Problem Auxiliary variables Results
TimeDurationOff__ThermalStandard
HydroEnergyOutput__HydroDispatch
TimeDurationOn__ThermalStandard
PowerSimulations Problem Expressions Results
ProductionCostExpression__RenewableDispatch
ProductionCostExpression__ThermalStandard
ActivePowerBalance__System
ProductionCostExpression__HydroDispatch
FuelConsumptionExpression__ThermalStandard
PowerSimulations Problem Parameters Results
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
PowerSimulations Problem Variables Results
StopVariable__ThermalStandard
ActivePowerReserveVariable__VariableReserve__ReserveUp__Reg_Up
ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R3
ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R1
ActivePowerVariable__ThermalStandard
ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R2
ActivePowerVariable__HydroDispatch
ActivePowerVariable__RenewableDispatch
OnVariable__ThermalStandard
StartVariable__ThermalStandard
ActivePowerReserveVariable__VariableReserve__ReserveDown__Reg_Down

Optimizer Stats

The optimizer summary is included

get_optimizer_stats(res)
1×21 DataFrame
Rowdetailed_statsobjective_valuetermination_statusprimal_statusdual_statussolver_solve_timeresult_counthas_valueshas_dualsobjective_boundrelative_gapdual_objective_valuesolve_timebarrier_iterationssimplex_iterationsnode_counttimed_solve_timetimed_calculate_aux_variablestimed_calculate_dual_variablessolve_bytes_allocsec_in_gc
BoolFloat64Int64Int64Int64Float64Int64BoolBoolMissingMissingMissingFloat64MissingMissingMissingFloat64Float64Float64Float64Float64
1false2.3571e6110NaN1falsefalsemissingmissingmissing0.578015missingmissingmissing0.6447110.4447320.0004674132.47632e70.0

Objective Function Value

get_objective_value(res)
2.357099085649291e6

Variable, Parameter, Auxillary Variable, Dual, and Expression Values

The solution value data frames for variables, parameters, auxillary variables, duals and expressions can be accessed using the read_ methods:

read_variables(res)
Dict{String, DataFrames.DataFrame} with 11 entries:
  "ActivePowerReserveVaria… => 24×52 DataFrame…
  "StopVariable__ThermalSt… => 24×55 DataFrame…
  "ActivePowerReserveVaria… => 24×52 DataFrame…
  "OnVariable__ThermalStan… => 24×55 DataFrame…
  "ActivePowerVariable__Hy… => 24×2 DataFrame…
  "ActivePowerReserveVaria… => 24×19 DataFrame…
  "StartVariable__ThermalS… => 24×55 DataFrame…
  "ActivePowerVariable__Th… => 24×55 DataFrame…
  "ActivePowerVariable__Re… => 24×30 DataFrame…
  "ActivePowerReserveVaria… => 24×18 DataFrame…
  "ActivePowerReserveVaria… => 24×17 DataFrame

Or, you can read a single parameter values for parameters that exist in the results.

list_parameter_names(res)
read_parameter(res, "ActivePowerTimeSeriesParameter__RenewableDispatch")
24×30 DataFrame
RowDateTime122_WIND_1324_PV_3312_PV_1102_PV_1101_PV_1324_PV_2313_PV_2104_PV_1101_PV_2309_WIND_1310_PV_2113_PV_1317_WIND_1314_PV_1324_PV_1103_PV_1303_WIND_1314_PV_2102_PV_2314_PV_3320_PV_1101_PV_3319_PV_1314_PV_4310_PV_1215_PV_1313_PV_1101_PV_4119_PV_1
DateTimeFloat64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64
12020-01-01T00:00:00713.20.00.00.00.00.00.00.00.0142.80.00.0795.10.00.00.0480.80.00.00.00.00.00.00.00.00.00.00.00.0
22020-01-01T01:00:00712.80.00.00.00.00.00.00.00.0139.10.00.0794.40.00.00.0634.90.00.00.00.00.00.00.00.00.00.00.00.0
32020-01-01T02:00:00708.40.00.00.00.00.00.00.00.0145.30.00.0773.60.00.00.0487.30.00.00.00.00.00.00.00.00.00.00.00.0
42020-01-01T03:00:00710.70.00.00.00.00.00.00.00.0144.80.00.0767.30.00.00.0432.70.00.00.00.00.00.00.00.00.00.00.00.0
52020-01-01T04:00:00701.40.00.00.00.00.00.00.00.0137.10.00.0752.20.00.00.0407.90.00.00.00.00.00.00.00.00.00.00.00.0
62020-01-01T05:00:00682.50.00.00.00.00.00.00.00.098.60.00.0719.40.00.00.0440.20.00.00.00.00.00.00.00.00.00.00.00.0
72020-01-01T06:00:00614.70.00.00.00.00.00.00.00.062.20.00.0655.30.00.00.0377.30.00.00.00.00.00.00.00.00.00.00.00.0
82020-01-01T07:00:00517.736.452.230.629.636.462.229.830.447.338.299.4594.647.444.829.8199.346.630.063.847.830.6179.246.438.2127.662.228.827.2
92020-01-01T08:00:00426.663.497.436.834.863.496.035.037.048.961.4126.2579.161.265.258.4110.667.036.4100.061.437.4248.067.061.4164.695.635.456.6
102020-01-01T09:00:00274.271.2118.038.036.471.4117.636.638.230.766.8133.2466.869.068.076.63.669.837.6121.869.838.8273.469.866.8172.2116.835.874.8
112020-01-01T10:00:0093.072.2132.037.836.472.2130.836.638.027.470.2134.2301.470.270.289.82.472.237.4129.070.438.6277.872.270.4173.2128.437.488.2
122020-01-01T11:00:006.370.6135.237.035.070.6134.235.437.260.969.8135.2110.767.070.298.056.272.436.6128.667.237.6263.272.469.6170.0133.238.691.6
132020-01-01T12:00:003.867.4131.036.235.467.4129.635.636.420.467.8133.078.967.268.294.091.570.235.8128.267.437.0267.470.267.6169.0128.837.890.6
142020-01-01T13:00:001.165.2125.435.235.065.2123.635.435.41.667.2126.6107.967.065.881.4103.067.634.8119.667.035.8258.867.667.0161.8122.835.879.2
152020-01-01T14:00:000.060.2109.631.031.060.4108.031.431.20.063.2110.222.359.462.659.439.764.230.8107.059.431.6227.864.262.8143.2106.831.459.0
162020-01-01T15:00:000.042.069.020.219.242.065.620.620.22.645.062.824.642.444.825.487.746.220.065.843.420.6151.646.243.882.262.820.427.8
172020-01-01T16:00:000.90.00.00.00.00.00.00.00.037.90.00.010.80.00.00.092.30.00.00.00.00.00.00.00.00.00.00.00.0
182020-01-01T17:00:00276.30.00.00.00.00.00.00.00.046.90.00.0243.20.00.00.089.40.00.00.00.00.00.00.00.00.00.00.00.0
192020-01-01T18:00:00272.90.00.00.00.00.00.00.00.024.00.00.0375.20.00.00.090.40.00.00.00.00.00.00.00.00.00.00.00.0
202020-01-01T19:00:00345.60.00.00.00.00.00.00.00.024.00.00.0568.40.00.00.081.10.00.00.00.00.00.00.00.00.00.00.00.0
212020-01-01T20:00:00411.70.00.00.00.00.00.00.00.027.40.00.0636.10.00.00.0172.90.00.00.00.00.00.00.00.00.00.00.00.0
222020-01-01T21:00:00376.60.00.00.00.00.00.00.00.06.50.00.0719.20.00.00.0326.90.00.00.00.00.00.00.00.00.00.00.00.0
232020-01-01T22:00:00561.30.00.00.00.00.00.00.00.01.30.00.0734.90.00.00.0256.70.00.00.00.00.00.00.00.00.00.00.00.0
242020-01-01T23:00:00568.40.00.00.00.00.00.00.00.00.10.00.0729.10.00.00.0141.10.00.00.00.00.00.00.00.00.00.00.00.0

Plotting

Take a look at the plotting capabilities in PowerGraphics.jl