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

sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")

Define a problem specification with an ProblemTemplate

You can create an empty template with:

template_uc = ProblemTemplate()

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

serialize_optimization_model(problem, save_path)

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

Optimizer Stats

The optimizer summary is included

get_optimizer_stats(res)

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…
  "ActivePowerReserveVaria… => 24×52 DataFrame…
  "StopVariable__ThermalSt… => 24×55 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")

Plotting

Take a look at the plotting capabilities in PowerGraphics.jl