ACOPF with PowerSimulations.jl using PowerModels.jl

Originally Contributed by: Clayton Barrows

Introduction

PowerSimulations.jl supports non-linear AC optimal power flow through a deep integration with PowerModels.jl. This example shows a single multi-period optimization of economic dispatch with a full representation of AC optimal power flow. However, since we use a case where generators are subject to minimum operating points, we need to also execute a unit commitment problem to provide the ACOPF with a valid commitment pattern. This example uses a Simulation with two DecisionModels to execute the UC-ACOPF workflow for a single period.

using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using PowerSystemCaseBuilder
using Dates
sim_folder = mktempdir(cleanup=true)

"/tmp/jl_yYPa5k"

We'll just use a suitable System that contains valid AC power flow parameters

sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")
transform_single_time_series!(sys, 1, Hour(1))

Since we'll be doing non-linear optimization, we need a solver that supports non-linear problems. Ipopt is quite good. And, we'll need a separate solver that can handle integer variables. So, we'll use HiGHS for the UC problem.

using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)
using HiGHS # mip solver
uc_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])

Here, we want do define an economic dispatch (linear generation decisions) with an ACOPF network representation. So, starting with the network, we can select from almost any of the endpoints on this tree:

# print_tree(PowerSimulations.PM.AbstractPowerModel)

First, we can setup a template with a suitable ACOPF network formulation, and formulations that represent each of the relevant device categories

ed_template = ProblemTemplate()
set_device_model!(ed_template, ThermalStandard, ThermalStandardDispatch)
set_device_model!(ed_template, PowerLoad, StaticPowerLoad)
set_device_model!(ed_template, Line, StaticBranch)
set_device_model!(ed_template, TapTransformer, StaticBranch)
set_device_model!(ed_template, Transformer2W, StaticBranch)
set_device_model!(ed_template, HVDCLine, HVDCDispatch)
set_network_model!(ed_template, NetworkModel(ACPPowerModel, use_slacks=true))

We also need to setup a UC template with a simplified network representation

uc_template = ProblemTemplate(DCPPowerModel)
set_device_model!(uc_template, ThermalStandard, ThermalBasicUnitCommitment)
set_device_model!(uc_template, PowerLoad, StaticPowerLoad)
set_device_model!(uc_template, Line, StaticBranch)
set_device_model!(uc_template, TapTransformer, StaticBranch)
set_device_model!(uc_template, Transformer2W, StaticBranch)
set_device_model!(uc_template, HVDCLine, HVDCDispatch)
set_service_model!(uc_template, VariableReserve{ReserveUp}, RangeReserve)

Now we can build a simulation to solve the UC, pass the commitment pattern to the ACOPF and then solve the ACOPF.

models = SimulationModels(
    decision_models=[
        DecisionModel(uc_template, sys, name="UC", optimizer=uc_solver),
        DecisionModel(
            ed_template,
            sys,
            name="ACOPF",
            optimizer=solver,
            initialize_model=false,
        ),
    ],
)
sequence = SimulationSequence(
    models=models,
    feedforwards=Dict(
        "ACOPF" => [
            SemiContinuousFeedforward(
                component_type=ThermalStandard,
                source=OnVariable,
                affected_values=[ActivePowerVariable, ReactivePowerVariable],
            ),
        ],
    ),
    ini_cond_chronology=InterProblemChronology(),
)

Note that in the above feedforward definition, the OnVariable for the ThermalStandard components is affecting both the ActivePowerVariable and the ReactivePowerVariable. This is the connection that restricts the ACOPF to only represent active and reactive power injections from the units that are committed in the UC problem. It's not guaranteed that the UC result generates an AC feasible initial condition, so this problem selects a particular period (initial_time) where the conditions are suitable.

sim = Simulation(
    name="UC-ACOPF",
    steps=12,
    models=models,
    sequence=sequence,
    simulation_folder=sim_folder,
)

build!(sim)

PowerSimulations.BuildStatusModule.BuildStatus.BUILT = 0

And solve it ...

execute!(sim, enable_progress_bar=false)

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

PowerSimulations.RunStatusModule.RunStatus.SUCCESSFUL = 0

And extract some results

results = SimulationResults(sim)
ac_results = get_decision_problem_results(results, "ACOPF")

slack_keys = [
    k for k in list_variable_keys(ac_results) if PSI.get_entry_type(k) ∈ [SystemBalanceSlackDown, SystemBalanceSlackUp]
]
slack_vars = Dict([k => read_realized_variables(ac_results, k) for k in slack_keys])

LoadError: MethodError: no method matching length(::PowerSimulations.VariableKey{PowerSimulations.SystemBalanceSlackDown, PowerSystems.Bus})
Closest candidates are:
  length(!Matched::Union{Base.KeySet, Base.ValueIterator}) at /opt/hostedtoolcache/julia/1.7.3/x64/share/julia/base/abstractdict.jl:58
  length(!Matched::Union{DataStructures.OrderedRobinDict, DataStructures.RobinDict}) at ~/.julia/packages/DataStructures/59MD0/src/ordered_robin_dict.jl:86
  length(!Matched::Union{DataStructures.SortedDict, DataStructures.SortedMultiDict, DataStructures.SortedSet}) at ~/.julia/packages/DataStructures/59MD0/src/container_loops.jl:322
  ...
Stacktrace:
  [1] length(g::Base.Generator{PowerSimulations.VariableKey{PowerSimulations.SystemBalanceSlackDown, PowerSystems.Bus}, PowerSimulations.var"#312#313"})
    @ Base ./generator.jl:50
  [2] _similar_shape(itr::Base.Generator{PowerSimulations.VariableKey{PowerSimulations.SystemBalanceSlackDown, PowerSystems.Bus}, PowerSimulations.var"#312#313"}, #unused#::Base.HasLength)
    @ Base ./array.jl:600
  [3] collect(itr::Base.Generator{PowerSimulations.VariableKey{PowerSimulations.SystemBalanceSlackDown, PowerSystems.Bus}, PowerSimulations.var"#312#313"})
    @ Base ./array.jl:723
  [4] read_realized_variables(res::PowerSimulations.SimulationProblemResults{PowerSimulations.DecisionModelSimulationResults}, variables::PowerSimulations.VariableKey{PowerSimulations.SystemBalanceSlackDown, PowerSystems.Bus}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PowerSimulations ~/.julia/packages/PowerSimulations/HKwUO/src/simulation/simulation_problem_results.jl:268
  [5] read_realized_variables(res::PowerSimulations.SimulationProblemResults{PowerSimulations.DecisionModelSimulationResults}, variables::PowerSimulations.VariableKey{PowerSimulations.SystemBalanceSlackDown, PowerSystems.Bus})
    @ PowerSimulations ~/.julia/packages/PowerSimulations/HKwUO/src/simulation/simulation_problem_results.jl:268
  [6] (::Main.__FRANKLIN_1802956.var"#9#10")(k::PowerSimulations.VariableKey{PowerSimulations.SystemBalanceSlackDown, PowerSystems.Bus})
    @ Main.__FRANKLIN_1802956 ./none:0
  [7] iterate
    @ ./generator.jl:47 [inlined]
  [8] collect(itr::Base.Generator{Vector{PowerSimulations.VariableKey{T, PowerSystems.Bus} where T<:PowerSimulations.VariableType}, Main.__FRANKLIN_1802956.var"#9#10"})
    @ Base ./array.jl:724

Plot the slack values

plot_results(slack_vars);