PowerFlow Validation Notes

[luke-kiernan]

In the process of validating the results of our power flow against that of matpower and now PSS/E, we've discovered some oddities of each. I'm documenting them here: partially for developers, partially for users.

Matpower

One thing that causes differences: PV buses where all generators are marked as unavailable. Matpower handles them correctly as PQ buses during the iterative Newton-Raphson process. However, after it has solved for the voltage, there's a bug in how it calculates (or perhaps fails to calculate) Q at such buses in pfsoln.m. Credit goes to Rodrigo for tracking that down. If you're comparing against matpower, pre-process the input file and change PV buses where all generators are unavailable to PQ. EDIT: looks like this was fixed at some point between 7.1 (bug described above present) and 8.0 (no such issue).

A second issue is that matpower always splits the total line charging susceptance, assigning half to the "from" bus and half to the "to" bus. In PowerSystems, we do the same thing for plain Lines, but not for various kinds of Transformers. (This might change: check this issue.) In the caseformat file, change the b value on all transformer branches to 0.0 and re-run: if the results then agree, this is your culprit.

Another transformer-related difference: matpower adds the line charging then multiplies by the tap ratio, whereas PowerNetworkMatrices does the opposite. In equations, matpower has Y_11 = (y_s + im*b/2)*c^2, where we have Y_11 = y_s*c^2 + im*b/2.

At the time of writing, PowerFlows.jl does not support nonzero reference angle. ACTIVSg10k has a nonzero reference angle, so avoid that system, or patch it in locally.

Script for comparing outputs against matpower

function normalize_angles!(data::PowerFlowData)
    # normalize angles to be within -pi to pi
    data.bus_angles[:, 1] .%= 2 * π
    data.bus_angles[:, 1] .+= 2 * π
    data.bus_angles[:, 1] .%= 2 * π
    # now they are in [0, 2pi)
    data.bus_angles[:, 1] .-= 2 * pi .* (data.bus_angles[:, 1] .> π)
end

# run in matlab: 
# myOpt = mpoption(mpoption, 'pf.nr.max_it', 100, 'pf.tol', 1e-12, 'verbose', 3)
# runpf('input.m', myOpt, 'output', 'results')

# parse the matpower results into Sienna's PowerFlowData struct
matpower_results = System("results.m")
matpower_data = PowerFlowData(ACPowerFlow(), matpower_results)

# run a powerflow on the same input
sys = System("input.m")
pf = ACPowerFlow{NewtonRaphsonACPowerFlow}()
sienna_data = PowerFlowData(pf, sys)
TOLERANCE = 1e-12 # may need to adjust: if 1e-12 fails, try 1e-11.
solve_powerflow!(sienna_data; maxIterations = 100, tol = TOLERANCE)

# some sanity checks
@assert matpower_data.bus_type[:, 1] == sienna_data.bus_type[:, 1]
PV_mask = matpower_data.bus_type[:, 1] .== (ACBusTypes.PV,)
@assert matpower_data.bus_magnitude[PV_mask, 1] == sienna_data.bus_magnitude[PV_mask, 1]

# compare
normalize_angles!(matpower_data)
normalize_angles!(sienna_data)
println("Sizes of differences: ")
println("\tVm: ", norm(matpower_data.bus_magnitude .- sienna_data.bus_magnitude, Inf))
println("\tVa: ", norm(matpower_data.bus_angles .- sienna_data.bus_angles, Inf))
println("\tGenerator p: ",
    norm(
        matpower_data.bus_activepower_injection .-
        sienna_data.bus_activepower_injection,
        Inf,
    ),
)
println("\tGenerator q: ",
  norm(
    matpower_data.bus_reactivepower_injection .-
    sienna_data.bus_reactivepower_injection,
    Inf,
  ),
)

[luke-kiernan]

PSSe

After running a power flow on PSS/e, I'll find that buses that were PV no longer match their voltage set point. That leads me to believe it's doing PV -> PQ switching on buses that hit their reactive power limits, even though we're passing -1 as the 7th option to FNSL (meaning: "ignore var limits"). Marck tried relaxing Q_min and Q_max, but it just hits the new reactive power bounds as well. [I do not have access to PSSe myself.]