code: def d2Sbr_dV2(self, Cbr, Ybr, V, lam):
          
          nb = len(V)
          nl = len(lam)
          ib = range(nb)
          il = range(nl)
  
          diaglam = csr_matrix((lam, (il, il)))
          diagV = csr_matrix((V, (ib, ib)))
  
          A = Ybr.H * diaglam * Cbr
          B = conj(diagV) * A * diagV
          D = csr_matrix( ((A * V) * conj(V), (ib, ib)) )
          E = csr_matrix( ((A.T * conj(V) * V), (ib, ib)) )
          F = B + B.T
          G = csr_matrix((ones(nb) / abs(V), (ib, ib)))
  
          Haa = F - D - E
          Hva = 1j * G * (B - B.T - D + E)
          Hav = Hva.T
          Hvv = G * F * G
  
          return Haa, Hav, Hva, Hvv
docs: Based on d2Sbr_dV2.m from MATPOWER by Ray Zimmerman, developed
          at PSERC Cornell. See U{http://www.pserc.cornell.edu/matpower/} for
          more information.
  
          @rtype: tuple
          @return: The 2nd derivatives of complex power flow w.r.t. voltage.