Returns the basin constraints for piece-wise linear gen cost
        variables.  CCV cost formulation expressed as Ay * x <= by.

        Based on makeAy.m from MATPOWER by C. E. Murillo-Sanchez, developed at
        PSERC Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more
        information.

def _pwl_gen_costs(self, generators, base_mva):
        """ Returns the basin constraints for piece-wise linear gen cost
        variables.  CCV cost formulation expressed as Ay * x <= by.

        Based on makeAy.m from MATPOWER by C. E. Murillo-Sanchez, developed at
        PSERC Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more
        information.
        """
        ng = len(generators)
        gpwl = [g for g in generators if g.pcost_model == PW_LINEAR]
#        nq = len([g for g in gpwl if g.qcost_model is not None])

        if self.dc:
            pgbas = 0        # starting index within x for active sources
            nq = 0           # number of Qg vars
#            qgbas = None     # index of 1st Qg column in Ay
            ybas = ng        # starting index within x for y variables
        else:
            pgbas = 0
            nq = ng
#            qgbas = ng + 1 # index of 1st Qg column in Ay
            ybas = ng + nq

        # Number of extra y variables.
        ny = len(gpwl)

        if ny == 0:
            return None, None

        # Total number of cost points.
        nc = len([co for gn in gpwl for co in gn.p_cost])
#        Ay = lil_matrix((nc - ny, ybas + ny))
        # Fill rows and then transpose.
        Ay = lil_matrix((ybas + ny, nc - ny))
        by = array([])

        j = 0
        k = 0
        for i, g in enumerate(gpwl):
            # Number of cost points: segments = ns-1
            ns = len(g.p_cost)

            p = array([x / base_mva for x, c in g.p_cost])
            c = array([c for x, c in g.p_cost])
            m = diff(c) / diff(p)        # Slopes for Pg (or Qg).

            if 0.0 in diff(p):
                raise ValueError, "Bad Pcost data: %s (%s)" % (p, g.name)
                logger.error("Bad Pcost data: %s" % p)

            b = m * p[:ns-1] - c[:ns-1] # rhs
            by = r_[by, b.T]

#            if i > ng:
#                sidx = qgbas + (i-ng) - 1       # this was for a q cost
#            else:
#                sidx = pgbas + i - 1            # this was for a p cost

            Ay[pgbas + i, k:k + ns - 1] = m

            # FIXME: Repeat for Q costs.

            # Now fill the y rows with -1's
            Ay[ybas + j, k:k + ns - 1] = -ones(ns-1)

            k += (ns - 1)
            j += 1

        y = Variable("y", ny)

        # Transpose Ay since lil_matrix stores in rows.
        if self.dc:
            ycon = LinearConstraint("ycon", Ay.T, None, by, ["Pg", "y"])
        else:
            ycon = LinearConstraint("ycon", Ay.T, None, by, ["Pg", "Qg","y"])

        return y, ycon