code: def d2Sbus_dV2(self, Ybus, V, lam):
          
          nb = len(V)
          ib = range(nb)
          Ibus = Ybus * V
          diaglam = csr_matrix((lam, (ib, ib)))
          diagV = csr_matrix((V, (ib, ib)))
  
          A = csr_matrix((lam * V, (ib, ib)))
          B = Ybus * diagV
          C = A * conj(B)
          D = Ybus.H * diagV
          E = diagV.conj() * (D * diaglam - csr_matrix((D * lam, (ib, ib))))
          F = C - A * csr_matrix((conj(Ibus), (ib, ib)))
          G = csr_matrix((ones(nb) / abs(V), (ib, ib)))
  
          Gaa = E + F
          Gva = 1j * G * (E - F)
          Gav = Gva.T
          Gvv = G * (C + C.T) * G
  
          return Gaa, Gav, Gva, Gvv
docs: Based on d2Sbus_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 power injection w.r.t. voltage.