code: def d2Sbus_dV2(Ybus, V, lam):
      
      n = len(V)
      Ibus = Ybus * V
      diaglam = spdiag(lam)
      diagV = spdiag(V)
  
      A = spmatrix(mul(lam, V), range(n), range(n))
      B = Ybus * diagV
      C = A * conj(B)
      D = Ybus.H * diagV
      E = conj(diagV) * (D * diaglam - spmatrix(D*lam, range(n), range(n)))
      F = C - A * spmatrix(conj(Ibus), range(n), range(n))
      G = spmatrix(div(matrix(1.0, (n, 1)), abs(V)), range(n), range(n))
  
      Gaa = E + F
      Gva = 1j * G * (E - F)
      Gav = Gva.T
      Gvv = G * (C + C.T) * G
  
      return Gaa, Gav, Gva, Gvv
docs: Computes 2nd derivatives of power injection w.r.t. voltage.