code: def _const_pf_constraints(self, gn, base_mva):
          ivl = array([i for i, g in enumerate(gn)
                       if g.is_load and (g.q_min != 0.0 or g.q_max != 0.0)])
          vl = [gn[i] for i in ivl]
          nvl = len(vl)
          ng = len(gn)
          Pg = array([g.p for g in vl]) / base_mva
          Qg = array([g.q for g in vl]) / base_mva
          Pmin = array([g.p_min for g in vl]) / base_mva
          Qmin = array([g.q_min for g in vl]) / base_mva
          Qmax = array([g.q_max for g in vl]) / base_mva
          for g in vl:
              if g.qmin != 0.0 and g.q_max != 0.0:
                  logger.error("Either Qmin or Qmax must be equal to zero for "
                  "each dispatchable load.")
          Qlim = (Qmin == 0.0) * Qmax + (Qmax == 0.0) * Qmin
          if any( abs(Qg - Pg * Qlim / Pmin) > 1e-6 ):
              logger.error("For a dispatchable load, PG and QG must be "
                           "consistent with the power factor defined by "
                           "PMIN and the Q limits.")
          if nvl > 0:
              xx = Pmin
              yy = Qlim
              pftheta = arctan2(yy, xx)
              pc = sin(pftheta)
              qc = -cos(pftheta)
              ii = array([range(nvl), range(nvl)])
              jj = r_[ivl, ivl + ng]
              Avl = csr_matrix(r_[pc, qc], (ii, jj), (nvl, 2 * ng))
              lvl = zeros(nvl)
              uvl = lvl
          else:
              Avl = zeros((0, 2 * ng))
              lvl = array([])
              uvl = array([])
          return LinearConstraint("vl", Avl, lvl, uvl, ["Pg", "Qg"])