Singular matrix error when solving simple 6 bus system

[rwl]

As reported by Dominik on the mailing list:

When including more than one generator in my grid, I get an error of creating a singular matrix:

[...]
  File "pypower\runpf.py", line 183, in runpf
    V, success, _ = newtonpf(Ybus, Sbus, V0, ref, pv, pq, ppopt)
  File "pypower\newtonpf.py", line 101, in newtonpf
    dx = -1 * spsolve(J, F)
  File "C:\WinPython-64bit-2.7.9.1\python-2.7.9.amd64\lib\site-packages\scipy\sparse\linalg\dsolve\linsolve.py", line 143, in spsolve
    b, flag, options=options)
RuntimeError: superlu failure (singular matrix?) at line 100 in file scipy\sparse\linalg\dsolve\SuperLU\SRC\dsnode_bmod.c

This does however not happen when solving the exact same case using MATPOWER after exporting it using 'savecase'.

def Testcase36():

    ## 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.0006, 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],
        [5, 1, 0, 0, 0, 0, 1, 1, -150, 0.4, 1, 1.1, 0.9],
        [6, 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],
        [5, 0.00104344028, -0, -0, -0, 0, 10, 1, 0.00104344028, 0.00104344028, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [6, 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],
        [3, 5, 2.253125, 0.879645943, 0, 0.190525589, 0.190525589, 0.190525589, 0, 0, 1, -360, 360],
        [5, 6, 1.3220625, 0.174947441, 0, 0.0997661265, 0.0997661265, 0.0997661265, 0, 0, 1, -360, 360],
    ])

    return ppc

[DMRWTH]

Just to give a short update, the issue could be solved by simply modelling all other generators as negative loads. This might pose problems if more complex OPF calculations are to be made, which is however not the case for me.

[Heiner92]

This is one major problem of PyPower. The Algorithm does not detect branches without anything connected.
What I'm trying to say is that when the last Bus of a strand has a total value of (0+j0)MVA the Jacobian cannot be inverted.
In the case above it would be Bus-4...

---

Simple Fix is to change all 0 values for P in ppc bus array to 1e-10

---

What should be done though is to automatically detect these branches/buses/generators and deactivate them. Before building the Jacobian.