Running Power Flow In The Loop with Unit Commitment

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

In this tutorial, you'll configure a unit commitment (UC) simulation that automatically calls an AC power flow solver from PowerFlows.jl at every dispatch interval.

You'll validate the committed dispatch against the full AC network model at each hour, check for branch overloads, and export the results to PSS/e format if needed for further analysis.

This tutorial builds on the tutorial for Running a Multi-Stage Production Cost Simulation.

Setup

Load the needed packages and define basic inputs:

using PowerSystemCaseBuilder
using PowerSimulations
using HydroPowerSimulations
using PowerFlows
using PowerSystems
using DataFrames
using HiGHS
using Dates
using Logging
const SCENARIO_NAME = "uc_with_pf"
const INLOOP_SYSTEM = "modified_RTS_GMLC_DA_sys"
const INLOOP_SIM_STEPS = 1

run_dir = joinpath(".", "uc_power_flow_in_the_loop_results")
mkpath(run_dir)
export_dir = joinpath(run_dir, "psse_exports")
mkpath(export_dir)

"./uc_power_flow_in_the_loop_results/psse_exports"

Build a test PowerSystems.System via PowerSystemCaseBuilder.build_system. We use modified_RTS_GMLC_DA_sys (24+ DA forecast steps):

sys = build_system(
    PSISystems,
    INLOOP_SYSTEM;
    skip_serialization = true,
    runchecks = false,
)

Configuring the Power Flow Solver

Create an PowerFlows.ACPowerFlow solver. We attach a PowerFlows.PSSEExportPowerFlow exporter so that we can automatically write a PSS/e .raw file for each solved interval.

psse_export = PSSEExportPowerFlow(;
    psse_version = :v33,
    export_dir = export_dir,
    overwrite = true,
)

power_flow_model = ACPowerFlow(; exporter = psse_export)

PowerFlows.ACPolarPowerFlow{PowerFlows.NewtonRaphsonACPowerFlow}(false, PowerFlows.PSSEExportPowerFlow(:v33, "./uc_power_flow_in_the_loop_results/psse_exports", "export", false, true), false, false, false, nothing, true, false, false, false, PowerNetworkMatrices.NetworkReduction[], 1, String[], false, false, Dict{Symbol, Any}())

Alternative: Fast Decoupled (FDNR) solver

PowerFlows 0.22 adds a fast-decoupled AC solver, selected via the solver type parameter of PowerFlows.ACPowerFlow. It plugs into power_flow_evaluation exactly like the default Newton-Raphson solver (same auxiliary variables and PSS/e exports) and matches its converged solution, but is often faster on large networks. Optionally hand off to an exact Newton solver for the final refinement:

using PowerFlows: FastDecoupledACPowerFlow, NewtonRaphsonACPowerFlow

fd_power_flow_model = ACPowerFlow{FastDecoupledACPowerFlow}(;
    exporter = psse_export,
    solver_settings = Dict{Symbol, Any}(
        :handoff_solver => NewtonRaphsonACPowerFlow,  # refine FD result with Newton-Raphson
        :handoff_tol => 1e-3,                         # FD-stage exit tolerance before handoff
    ),
)

Use it in place of power_flow_model below. Omit solver_settings for pure fast decoupled, or use the FastDecoupledFixed alias for the formulation-agnostic frozen-Jacobian variant.

Building the UC Problem Template

Create a ProblemTemplate with a NetworkModel that uses PTDFPowerModel and power_flow_evaluation set to the solver we just configured. Assign device formulations with set_device_model!. This is the key step that enables power flow in the loop:

template_uc = ProblemTemplate(
    NetworkModel(
        PTDFPowerModel;
        use_slacks = true,
        power_flow_evaluation = power_flow_model,
    ),
)

set_device_model!(
    template_uc,
    ThermalStandard,
    ThermalStandardUnitCommitment,
)
set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch)
set_device_model!(template_uc, RenewableNonDispatch, FixedOutput)
set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
set_device_model!(template_uc, HydroDispatch, HydroDispatchRunOfRiver)
set_device_model!(template_uc, Line, StaticBranchUnbounded)
set_device_model!(template_uc, Transformer2W, StaticBranchUnbounded)
set_device_model!(template_uc, MonitoredLine, StaticBranch)

use_slacks = true allows the simulation to remain feasible when there is a small mismatch between the PTDF-based UC network model and the full AC power flow.

Building and Executing the Simulation

Follow the same pattern as in the tutorial for Running a Multi-Stage Production Cost Simulation: package the UC DecisionModel in SimulationModels, define a SimulationSequence with InterProblemChronology, construct a Simulation, then build! and execute! it for one simulation step. The in-step horizon on modified_RTS_GMLC_DA_sys provides 24 hourly realized rows.

solver = optimizer_with_attributes(
    HiGHS.Optimizer,
    "log_to_console" => false,
    "mip_rel_gap" => 0.05,
    "time_limit" => 900.0,
)

models = SimulationModels(;
    decision_models = [
        DecisionModel(
            template_uc,
            sys;
            name = "UC",
            optimizer = solver,
            store_variable_names = true,
        ),
    ],
)

sequence = SimulationSequence(;
    models = models,
    ini_cond_chronology = InterProblemChronology(),
)

sim = Simulation(;
    name = SCENARIO_NAME,
    steps = INLOOP_SIM_STEPS,
    models = models,
    sequence = sequence,
    simulation_folder = run_dir,
)

build!(sim; console_level = Logging.Error, file_level = Logging.Error)
execute!(sim; enable_progress_bar = true)

InfrastructureSystems.Simulation.RunStatusModule.RunStatus.SUCCESSFULLY_FINALIZED = 0

Loading Simulation Results

Load the simulation results with SimulationResults and extract the UC problem results with get_decision_problem_results:

sim_results = SimulationResults(sim)
uc_results = get_decision_problem_results(sim_results, "UC")

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

End: 2020-01-01T00:00:00

Resolution: 60 minutes

Notice the "UC Problem Auxiliary variables Results" table, which lists the active and reactive power flow and bus voltage magnitude and angle results from the AC power flow (e.g., PowerFlowBranchReactivePowerFromTo__Line, PowerFlowVoltageMagnitude__ACBus). These are not output when a UC problem is run alone.

Power Flow Auxiliary Variable Types

The AC power flow writes its solved quantities into auxiliary variables, all subtypes of PowerFlowAuxVariableType. They are grouped by what each is indexed on:

PowerFlowAuxVariableType
├─ BranchFlowAuxVariableType              # per-AC-branch flows (natural-units powers)
│    ├─ PowerFlowBranchActivePowerFromTo / …ToFrom / …Loss
│    └─ PowerFlowBranchReactivePowerFromTo / …ToFrom
├─ PowerFlowHVDCAuxVariableType           # per-HVDC-component, from PowerFlows.get_hvdc_results
│    ├─ PowerFlowHVDCActivePower/ReactivePower {FromTo,ToFrom} + PowerFlowHVDCActivePowerLoss
│    ├─ PowerFlowHVDCDCCurrent / DCVoltageFrom / DCVoltageTo
│    ├─ PowerFlowLCC{Rectifier,Inverter}Tap / RectifierDelayAngle / InverterExtinctionAngle
│    └─ PowerFlowConverterDCPower / ReactivePower / DCVoltage
├─ bus-indexed:   PowerFlowVoltageAngle, PowerFlowVoltageMagnitude, PowerFlowHVDCNetPower,
│                 PowerFlowLossFactors, PowerFlowVoltageStabilityFactors
└─ control-device-indexed:  PowerFlowTapRatio (TapTransformer),
                 PowerFlowSwitchedShuntSusceptance (SwitchedAdmittance),
                 PowerFlowFACTSReactivePower (FACTSControlDevice)

PowerFlowHVDCNetPower is a per-bus net-injection quantity, which is why it is a direct subtype of PowerFlowAuxVariableType and not a PowerFlowHVDCAuxVariableType (those are all per-component).

PTDF UC Flows vs. AC Power Flow In the Loop

Now, we'll compare the PTDF UC flows to the AC power flow results.

The UC stage optimizes flows using the PTDF network model, and the line flows are not variables; they are recorded in the PTDFBranchFlow__* expressions.

First, load in the PTDF branch flows for one line with read_realized_expression:

ptdf_flows = read_realized_expression(uc_results, "PTDFBranchFlow__Line")
example_line = first(unique(ptdf_flows.name))
ptdf_line = filter(row -> row.name == example_line, ptdf_flows)

After each in-step hour, the AC power flow writes branch flows into auxiliary variables such as PowerFlowBranchActivePowerFromTo__Line. Load those flows with read_realized_aux_variable and compare PTDF expression flows to AC flows for our selected line:

pf_flows_ft = read_realized_aux_variable(
    uc_results,
    "PowerFlowBranchActivePowerFromTo__Line",
)
pf_line = filter(row -> row.name == example_line, pf_flows_ft)

innerjoin(
    rename(select(ptdf_line, :DateTime, :value), :value => :PTDF_MW),
    rename(select(pf_line, :DateTime, :value), :value => :AC_PF_MW);
    on = :DateTime,
)

With use_slacks = true in the template, small differences between PTDF and AC flows are expected — the UC model can absorb minor mismatches. The AC auxiliary flows should track the PTDF values closely when the network is well conditioned.

Checking for Branch Overloads

Now, we will run a post-hoc check on our UC results to scan every interval for branches whose AC flow exceeds the thermal rating.

First, build a helper to extract branch flow limits from the PowerSystems.System — name in the auxiliary results matches the branch name on each PowerSystems.ACBranch. Limits are scaled with PowerSystems.get_base_power. AC lines and transformers use PowerSystems.get_rating; HVDC branches use active-power limits instead.

function _branch_flow_limit_mw(branch, sys)
    base_power = get_base_power(sys)
    try
        if branch isa TwoTerminalVSCLine
            return get_rating(branch) * base_power
        elseif branch isa TwoTerminalHVDC
            return max(
                get_active_power_limits_from(branch).max,
                get_active_power_limits_to(branch).max,
            )
        elseif branch isa GenericArcImpedance
            return get_max_flow(branch) * base_power
        else
            return get_rating(branch) * base_power
        end
    catch e
        e isa MethodError || rethrow(e)
        @warn "Could not get rating for $(typeof(branch)): $e — treating as unlimited"
        return Inf
    end
end

_branch_flow_limit_mw (generic function with 1 method)

Next, build a lookup for each PowerSystems.ACBranch keyed by branch name using PowerSystems.get_components and get_name:

ratings = Dict{String, Float64}()
for b in get_components(ACBranch, sys)
    ratings[get_name(b)] = _branch_flow_limit_mw(b, sys)
end

Then, scan every interval for branches whose AC flow exceeds the thermal rating:

overloads = DataFrame(;
    DateTime = Dates.DateTime[],
    name = String[],
    P_from_to = Float64[],
    rating = Float64[],
)

for row in eachrow(pf_flows_ft)
    rating = get(ratings, row.name, Inf)
    if abs(row.value) > rating
        push!(overloads, (row.DateTime, row.name, row.value, rating))
    end
end

overloads

Notice we have multiple overloads in this example. Each row identifies a congestion event: the interval (DateTime), branch name (name), actual AC flow in MW (P_from_to), and thermal rating (rating). In-loop AC power flow records realized flows but does not add thermal-limit constraints to the UC optimization. The UC model's network representation could be too loose if overloads appear — add transmission constraints, re-run UC, and repeat the AC power flow check until the table is empty.

PSS/e Export Files

The PowerFlows.PSSEExportPowerFlow exporter writes under export_dir (here, results/psse_exports). After running, we have one folder per solved interval in the 24-hour in-step horizon, plus the 24-hour lookahead:

psse_export_folders = readdir(export_dir)
length(psse_export_folders)

48

Each folder contains the PSS/e .raw file for that interval, available for further analysis:

first_subdir = first(sort(psse_export_folders))
readdir(joinpath(export_dir, first_subdir))

2-element Vector{String}:
 "export_1_1.raw"
 "export_1_1_export_metadata.json"

Next Steps

• Tighten the UC network model with transmission constraints if needed until overloads are removed — see the formulation library Network and PowerSystems.Branch Formulations pages.
• See the Solving a Power Flow tutorial in PowerFlows.jl for standalone power flow examples.