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"])