code: def dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St):
      
      dAf_dPf = spdiag(2 * Sf.real())
      dAf_dQf = spdiag(2 * Sf.imag())
      dAt_dPt = spdiag(2 * St.real())
      dAt_dQt = spdiag(2 * St.imag())
  
      # Partial derivative of apparent power magnitude w.r.t voltage
      # phase angle.
      dAf_dVa = dAf_dPf * dSf_dVa.real() + dAf_dQf * dSf_dVa.imag()
      dAt_dVa = dAt_dPt * dSt_dVa.real() + dAt_dQt * dSt_dVa.imag()
      # Partial derivative of apparent power magnitude w.r.t. voltage
      # amplitude.
      dAf_dVm = dAf_dPf * dSf_dVm.real() + dAf_dQf * dSf_dVm.imag()
      dAt_dVm = dAt_dPt * dSt_dVm.real() + dAt_dQt * dSt_dVm.imag()
  
      return dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm
docs: Partial derivatives of squared flow magnitudes w.r.t voltage.
  
          Computes partial derivatives of apparent power w.r.t active and
          reactive power flows.  Partial derivative must equal 1 for lines
          with zero flow to avoid division by zero errors (1 comes from
          L'Hopital).