code: def d2Sbr_dV2(Cbr, Ybr, V, lam):
      
      nb = len(V)
  
      diaglam = spdiag(lam)
      diagV = spdiag(V)
  
      A = Ybr.H * diaglam * Cbr
      B = conj(diagV) * A * diagV
      D = spdiag(mul((A*V), conj(V)))
      E = spdiag(mul((A.T * conj(V)), V))
      F = B + B.T
      G = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))
  
      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: Computes 2nd derivatives of complex power flow w.r.t. voltage.