Solving a Power Flow

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

In this tutorial, you'll solve power flows on a 5-bus test system using three different solvers and compare their results.

Building a System

To get started, load the needed packages. We're using a standard test system and want to keep output clean, so we adjust the logging settings to filter out a few precautionary warnings.

using PowerSystemCaseBuilder
using PowerFlows
using PowerSystems
using Logging

disable_logging(Logging.Warn)
LogLevel(1001)

Create a System from PowerSystemCaseBuilder.jl:

sys = build_system(MatpowerTestSystems, "matpower_case5_sys")
System
Property Value
Name
Description
System Units Base SYSTEM_BASE
Base Power 100.0
Base Frequency 60.0
Num Components 28
Static Components
Type Count
ACBus 5
Arc 6
Area 1
Line 5
LoadZone 1
PhaseShiftingTransformer 2
PowerLoad 3
ThermalStandard 5

DC Power Flow

DCPowerFlow solves for bus voltage angles using the bus admittance matrix, then computes branch flows from the angle differences. Create a DCPowerFlow solver:

pf_dc = DCPowerFlow()
DCPowerFlow(nothing, PowerNetworkMatrices.NetworkReduction[], 1, String[], false, false)

Solve the power flow with solve_power_flow. For DC methods, pass a FlowReporting mode so flows are reported on the same arc basis as elsewhere in the package:

dc_results = solve_power_flow(pf_dc, sys, FlowReporting.ARC_FLOWS)
Dict{String, Dict{String, DataFrame}} with 1 entry:
  "1" => Dict("flow_results"=>6×10 DataFrame

The result is a Dict{String, Dict{String, DataFrame}}. The outer key is the time step name: "1". The inner dictionary stores the power flow results at that time step: "bus_results" for bus data and "flow_results" for AC line data. (There is also a third key, "lcc_results", for HVDC lines, but this system contains no such components, so the matching dataframe will be empty.) Inspect "bus_results":

dc_results["1"]["bus_results"]
bus_number Vm θ P_gen P_load P_net Q_gen Q_load Q_net
Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
1 1.0 0.05813712342707429 210.0 0.0 210.0 0.0 0.0 0.0
2 1.0 -0.010697105096729564 0.0 300.0 -300.0 0.0 0.0 0.0
3 1.0 -0.004752965105551588 324.498 300.0 24.497999999999998 0.0 0.0 0.0
4 1.0 0.0 0.0 400.0 -400.0 0.0 0.0 0.0
10 1.0 0.07261406474007273 470.694 0.0 470.694 0.0 0.0 0.0

Notice that Vm (voltage magnitude) is 1.0 for all buses, and Q_gen and Q_load are 0. This is expected for DC power flow, which assumes flat voltage magnitudes and ignores reactive power.

dc_results["1"]["flow_results"]
flow_name bus_from bus_to P_from_to Q_from_to P_to_from Q_to_from P_losses Q_losses angle_difference
String Int64 Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
1-2 1 2 244.9616680920614 0.0 -244.9616680920614 0.0 1.6861747492479113 0.0 0.06883422852380386
1-4 1 4 191.2405399235394 0.0 -191.2405399235394 0.0 1.1118175009515048 0.0 0.05813712342707429
1-10 1 10 -226.20220801560063 0.0 226.20220801560063 0.0 0.32747160903125155 0.0 -0.01447694131299844
2-3 2 3 -55.03833190793862 0.0 55.03833190793862 0.0 0.03271555417545088 0.0 -0.005944139991177976
3-4 3 4 -30.5403319079386 0.0 30.5403319079386 0.0 0.03054095074902095 0.0 -0.004752965105551588
4-10 4 10 -244.49179198439933 0.0 244.49179198439933 0.0 1.7753542195279608 0.0 -0.07261406474007273

Likewise, Q_from_to and Q_to_from (reactive power flow on the line) are zero, for all lines.

PTDF DC Power Flow

PTDFDCPowerFlow computes branch flows directly from bus power injections using the Power Transfer Distribution Factor matrix, without solving for voltage angles as an intermediate step. (This means we can omit the angle computation in contexts where we only care about line flows, though we don't have that option implemented here.) Create a PTDFDCPowerFlow solver:

pf_ptdf = PTDFDCPowerFlow()
PTDFDCPowerFlow(nothing, false, PowerNetworkMatrices.NetworkReduction[], 1, String[], false)

As before, solve the power flow with solve_power_flow:

ptdf_results = solve_power_flow(pf_ptdf, sys, FlowReporting.ARC_FLOWS)
Dict{String, Dict{String, DataFrame}} with 1 entry:
  "1" => Dict("flow_results"=>6×10 DataFrame

Look at the bus results:

ptdf_results["1"]["bus_results"]
bus_number Vm θ P_gen P_load P_net Q_gen Q_load Q_net
Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
1 1.0 0.05813712342707429 210.0 0.0 210.0 0.0 0.0 0.0
2 1.0 -0.010697105096729564 0.0 300.0 -300.0 0.0 0.0 0.0
3 1.0 -0.004752965105551588 324.498 300.0 24.497999999999998 0.0 0.0 0.0
4 1.0 0.0 0.0 400.0 -400.0 0.0 0.0 0.0
10 1.0 0.07261406474007273 470.694 0.0 470.694 0.0 0.0 0.0

The results match DCPowerFlow, as they should: the two are mathematically equivalent. For very large systems where forming the full PTDF matrix would be too expensive, consider vPTDFDCPowerFlow, which computes the same results without storing the dense matrix.

AC Power Flow

Create an ACPowerFlow solver:

pf_ac = ACPowerFlow()
ACPowerFlow{NewtonRaphsonACPowerFlow}(false, nothing, false, false, nothing, true, false, false, false, PowerNetworkMatrices.NetworkReduction[], 1, String[], false, Dict{Symbol, Any}())

Solve the power flow:

ac_results = solve_power_flow(pf_ac, sys)
Dict{String, DataFrame} with 3 entries:
  "flow_results" => 6×10 DataFrame…
  "bus_results"  => 5×9 DataFrame…
  "lcc_results"  => 0×13 DataFrame

AC results are returned as a flat Dict{String, DataFrame}, with the same keys as before: "bus_results", "flow_results" (AC lines), and "lcc_results" (HVDC lines). (We don't support multi-period AC power flows yet.) Look at the bus results:

ac_results["bus_results"]
bus_number Vm θ P_gen P_load P_net Q_gen Q_load Q_net
Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
1 1.07762 0.05377973212862824 210.0 0.0 210.0 155.86199106764073 0.0 155.86199106764073
2 1.0841946221315348 8.152349463538256e-5 0.0 300.0 -300.0 0.0 98.61 -98.61
3 1.1 0.006160589908746271 324.498 300.0 24.497999999999998 518.182149466319 98.61 419.57214946631893
4 1.06414 0.0 2.2015611506636112 400.0 -397.7984388493364 -116.00444489674346 131.47 -247.47444489674345
10 1.06907 0.06665597012564134 470.694 0.0 470.694 -164.31698920802899 0.0 -164.31698920802899

Notice that Vm now varies across buses (not all 1.0), and Q_gen has non-zero values.

Look at the line flows:

ac_results["flow_results"]
flow_name bus_from bus_to P_from_to Q_from_to P_to_from Q_to_from P_losses Q_losses angle_difference
String Int64 Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
1-2 1 2 219.04815651765563 -41.538343627190756 -217.84617277462783 52.726508142521936 1.2019837430278013 11.188164515331177 0.053698208633992855
1-4 1 4 206.03176235568074 32.25263490048205 -204.89263500836265 -21.615578151334624 1.139127347318114 10.637056749147423 0.05377973212862824
1-10 1 10 -215.0798358132005 165.14902875652393 215.48841730273134 -164.66348802499908 0.40858148953084594 0.4855407315248472 -0.012876237997013103
2-3 2 3 -82.153855250951 -151.33695655178607 82.4232748938078 151.82170300927572 0.2694196428567941 0.4847464574896376 -0.006079066414110889
3-4 3 4 -57.925274893807945 267.75044645704486 60.60738862294479 -241.67008678435033 2.682113729136848 26.08035967269453 0.006160589908746271
4-10 4 10 -253.51319246391864 15.812083986954987 255.20569169372328 0.3464988169711386 1.6924992298046337 16.158582803926123 -0.06665597012564134

Q_from_to and Q_to_from now show reactive power flows, and P_from_to differs from P_to_from due to losses.

When AC Power Flow Fails

Unlike DC power flow, AC power flow is iterative and not guaranteed to converge. Systems with high impedance lines, poor initial voltage profiles, or insufficient reactive power support can cause the solver to fail. When this happens, solve_power_flow returns missing: you'll also see a logged error. If you encounter convergence failures, consider using a more robust solver such as TrustRegionACPowerFlow or RobustHomotopyPowerFlow.