Returns a linear constraint enforcing constant power factor for
        dispatchable loads.

        The power factor is derived from the original value of Pmin and either
        Qmin (for inductive loads) or Qmax (for capacitive loads). If both Qmin
        and Qmax are zero, this implies a unity power factor without the need
        for an additional constraint.

def _const_pf_constraints(self, gn, base_mva):
        """ Returns a linear constraint enforcing constant power factor for
        dispatchable loads.

        The power factor is derived from the original value of Pmin and either
        Qmin (for inductive loads) or Qmax (for capacitive loads). If both Qmin
        and Qmax are zero, this implies a unity power factor without the need
        for an additional constraint.
        """
        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

        # At least one of the Q limits must be zero (corresponding to Pmax==0).
        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.")

        # Initial values of PG and QG must be consistent with specified power
        # factor. This is to prevent a user from unknowingly using a case file
        # which would have defined a different power factor constraint under a
        # previous version which used PG and QG to define the power factor.
        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.")

        # Make Avl, lvl, uvl, for lvl <= Avl * r_[Pg, Qg] <= uvl
        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"])