Losses incorrectly calculated

[rwl]

As reported by Dominik on the mailing list:

When running a Powerflow in a very simple network of only 4 busses PYPOWER delivers the exact same result as when running the case in MATPOWER, however in its output it estimates an enormous power loss in the first branch, which is obviously false. (0.79MW compared to a value below 1kW).

def Testcase99():

    ## PYPOWER Case Format : Version 2
    ppc = {'version': '2'}

    ##-----  Power Flow Data  -----##
    ## system MVA base
    ppc['baseMVA'] = 10

    ## bus data
    # bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
    ppc['bus'] = array([
        [1, 3, 0, 0, 0, 0, 1, 1.02, 0, 20, 1, 1.06, 0.98],
        [2, 1, 0, 0, 0, 0, 1, 1, -150, 0.4, 1, 1.1, 0.9],
        [3, 1, 0.000425005103, 0, 0, 0, 1, 1, -150, 0.4, 1, 1.1, 0.9],
        [4, 1, 0, 0, 0, 0, 1, 1, -150, 0.4, 1, 1.1, 0.9],
    ])

    ## generator data
    # bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
    ppc['gen'] = array([
        [1, 0, 0, 10, -10, 1.02, 10, 1, 10, -10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [2, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [3, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [4, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    ])

    ## branch data
    # fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
    ppc['branch'] = array([
        [1, 2, 2.52, 5.17949805, 0, 0.1, 0.1, 0.1, 1, 150, 1, -360, 360],
        [2, 3, 2.253125, 0.879645943, 0, 0.190525589, 0.190525589, 0.190525589, 0, 0, 1, -360, 360],
        [3, 4, 0.721125, 0.0954258769, 0, 0.0997661265, 0.0997661265, 0.0997661265, 0, 0, 1, -360, 360],
    ])

    return ppc

[DMRWTH]

I corrected this issue by updating algorithms used for loss calculation to a transcription of what is implemented in the latest version of MATPOWER. Area differentiation is however likely not working as it is of no concern for my work.

from sys import stdout

from scipy.sparse import bsr_matrix as bsr_matrix

from numpy import \
    ones, zeros, r_, sort, exp, pi, diff, arange, min, \
    argmin, argmax, logical_or, real, imag, any, array, diag, conj
from numpy import flatnonzero as find

from pypower.idx_bus import BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, \
    VM, VA, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN, REF
from pypower.idx_gen import GEN_BUS, PG, QG, QMAX, QMIN, GEN_STATUS, \
    PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN
from pypower.idx_brch import F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, \
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST

from pypower.isload import isload
from pypower.run_userfcn import run_userfcn
from pypower.ppoption import ppoption

def printpf(baseMVA, bus=None, gen=None, branch=None, f=None, success=None,
            et=None, fd=None, ppopt=None):
    [...]

    ## create map of external bus numbers to bus indices
    i2e = bus[:, BUS_I].astype(int)
    e2i = zeros(max(i2e) + 1, int)
    e2i[i2e] = arange(bus.shape[0])

    ## sizes of things
    nb = bus.shape[0]      ## number of buses
    nl = branch.shape[0]   ## number of branches
    ng = gen.shape[0]      ## number of generators

    ## zero out some data to make printout consistent for DC case
    if isDC:
        bus[:, r_[QD, BS]]          = zeros((nb, 2))
        gen[:, r_[QG, QMAX, QMIN]]  = zeros((ng, 3))
        branch[:, r_[BR_R, BR_B]]   = zeros((nl, 2))

    ## parameters
    ties = find(bus[e2i[branch[:, F_BUS].astype(int)], BUS_AREA] !=
                   bus[e2i[branch[:, T_BUS].astype(int)], BUS_AREA])

    tap = ones(nl)                           ## default tap ratio = 1 for lines
    xfmr = find(branch[:, TAP])           ## indices of transformers

    for nonzerotab in xfmr:                                                         
            tap[nonzerotab] = branch[nonzerotab, TAP]
    ##connection matrices       
    Cf = bsr_matrix((branch[:,BR_STATUS].astype(int),(array(range(nl)),array(e2i[branch[:,F_BUS].astype(int)]))))
    Cf = sparse.hstack([Cf,zeros((Cf.shape[0],1))])

    Ct = bsr_matrix((branch[:,BR_STATUS].astype(int),(array(range(nl)),array(e2i[branch[:,T_BUS].astype(int)]))))

    tap = tap * exp(1j * pi / 180 * branch[:, SHIFT]) ## add phase shifters
    nzld = find((bus[:, PD] != 0.0) | (bus[:, QD] != 0.0))
    sorted_areas = sort(bus[:, BUS_AREA])
    ## area numbers
    s_areas = sorted_areas[r_[1, find(diff(sorted_areas)) + 1]]
    nzsh = find((bus[:, GS] != 0.0) | (bus[:, BS] != 0.0))
    allg = find( ~isload(gen) )
    ong  = find( (gen[:, GEN_STATUS] > 0) & ~isload(gen) )
    onld = find( (gen[:, GEN_STATUS] > 0) &  isload(gen) )

    out = find(branch[:, BR_STATUS] == 0)        ## out-of-service branches
    nout = len(out)

    A = diag(1.0/tap)*Cf - Ct
    print A
    Ysc = 1.0/(branch[:,BR_R] - 1j * branch[:,BR_X])   
    V = bus[:, VM] * exp(1j * pi / 180 * bus[:, VA])
    Vtrans = V.reshape(-1,1)
    Vdrop = A*Vtrans
    loss = ones(nl)*(0+0j)

    if isDC:
        loss = zeros(nl)
    else:
        '''loss = baseMVA * abs(V[e2i[ branch[:, F_BUS].astype(int) ]] / tap -
                             V[e2i[ branch[:, T_BUS].astype(int) ]])**2 / \
                    (branch[:, BR_R] - 1j * branch[:, BR_X])'''
        for lineloss in range(nl):
            loss[lineloss] = (baseMVA*Ysc[lineloss]*Vdrop[lineloss]*conj(Vdrop[lineloss]))[0,0]

    fchg = abs(V[e2i[ branch[:, F_BUS].astype(int) ]] / tap)**2 * branch[:, BR_B] * baseMVA / 2
    tchg = abs(V[e2i[ branch[:, T_BUS].astype(int) ]]      )**2 * branch[:, BR_B] * baseMVA / 2
    loss[out] = zeros(nout)
    fchg[out] = zeros(nout)
    tchg[out] = zeros(nout)

    ##----- print the stuff -----
    [...]

[HeinerFrueh]

Hey there, I just wanted to let everybody who ran into similar problems know how to correct the false values for transformer losses:

the incorrect line in the printpf script is line 150: "tap = tap * exp(1j * pi / 180 * branch[:, SHIFT]) ## add phase shifters" it should say "tap = tap * exp( - 1j * pi / 180 * branch[:, SHIFT])" instead (note the negative angle!)

this only has an effect on transformers which are also phase shifting... I hope it will help somebody ;)

[Heiner92]

Are there still people around who maintain this project? Everytime I install PYPOWER on a new machine I have to make this change myself, which is starting to bother me. The fix is simple: just add the "-" ... I would very much appreciate it if this was fixed :)

the incorrect line in the printpf script is line 150: "tap = tap exp(1j pi / 180 branch[:, SHIFT]) ## add phase shifters" it should say "tap = tap exp( - 1j pi / 180 branch[:, SHIFT])" instead (note the negative angle!)