Inconsistent power flow results of function `pf_large`

[jinningwang]

Use this issue as a discussion board

[jinningwang]

Issue description: function pf_large give different results compared with ANDES

Solved voltage from pf_large and ANDES: case4gs.m: array([1. , 0.9805, 0.9666, 1.02 ]) array([1. , 0.9824, 0.969 , 1.02 ])

case5.m array([1. , 0.9892, 1. , 1. , 1. ]) array([1. , 0.9893, 1. , 1. , 1. ])

case9.m array([1.04 , 1.025 , 1.025 , 0.9957, 0.9662, 1.0041, 0.9792, 0.9974, 0.9507]) array([1.04 , 1.025 , 1.025 , 1.0258, 1.0127, 1.0324, 1.0159, 1.0258, 0.9956])

Function:

case = andes.get_case('matpower/case5.m')

ss = andes.load(case,
                default_config=True)

Y = matrix(ss.Line.build_y())

nbus = ss.Bus.n
P = np.zeros(nbus)
Q = np.zeros(nbus)
V = np.ones(nbus)
theta = np.zeros(nbus)

pq_loc = ss.Bus.idx2uid(ss.PQ.bus.v)
pv_loc = ss.Bus.idx2uid(ss.PV.bus.v)
slack_loc = ss.Bus.idx2uid(ss.Slack.bus.v)
gen_loc = np.concatenate([pv_loc, slack_loc])
non_gen_loc = np.setdiff1d(np.arange(ss.Bus.n), gen_loc)
non_slack_loc = np.setdiff1d(np.arange(ss.Bus.n), slack_loc)
shunt_loc = ss.Bus.idx2uid(ss.Shunt.bus.v)

def pf_large(x):

    # these are the arrays for each bus; 
    #   elements corresponding to bus numbers in ascending order
    P[:] = 0
    Q[:] = 0
    V[:] = 1
    theta[:] = 0

    # PQ
    P[pq_loc] -= ss.PQ.p0.v
    Q[pq_loc] -= ss.PQ.q0.v

    # PV
    V[pv_loc] = ss.PV.v0.v
    P[pv_loc] += ss.PV.p0.v

    # Slack
    theta[slack_loc] = ss.Slack.a0.v
    V[slack_loc] = ss.Slack.v0.v

    # retrieve unknowns from `x`
    # unknowns in `x` are grouped by theta and then V
    #   note that the `non_gen_loc` can have random order
    #   no sorting is needed, because x will be in the same order

    theta[non_slack_loc] = x[:len(non_slack_loc)]
    V[non_gen_loc] = x[len(non_slack_loc):]

    # shunt elements
    P[shunt_loc] += ss.Shunt.g.v * V[shunt_loc] ** 2
    Q[shunt_loc] += ss.Shunt.b.v * V[shunt_loc] ** 2

    Vc = V * np.exp(1j * theta)

    # calculate the power injection into the network
    #   note that power leaves the bus into the network
    S = np.diag(Vc) @ np.conj(Y @ Vc)

    # power leaving each bus via the lines 
    #   minus power injection at each bus shall equal to 0
    Pmis = np.real(S) - P
    Qmis = np.imag(S) - Q

    return np.concatenate([Pmis[non_slack_loc],
                           Qmis[non_gen_loc]])