Python numpy.find函数代码示例

def makePTDF(baseMVA, bus, branch, slack=None):
    """Builds the DC PTDF matrix for a given choice of slack.

    Returns the DC PTDF matrix for a given choice of slack. The matrix is
    C{nbr x nb}, where C{nbr} is the number of branches and C{nb} is the
    number of buses. The C{slack} can be a scalar (single slack bus) or an
    C{nb x 1} column vector of weights specifying the proportion of the
    slack taken up at each bus. If the C{slack} is not specified the
    reference bus is used by default.

    For convenience, C{slack} can also be an C{nb x nb} matrix, where each
    column specifies how the slack should be handled for injections
    at that bus.

    @see: L{makeLODF}

    @author: Ray Zimmerman (PSERC Cornell)
    """
    ## use reference bus for slack by default
    if slack is None:
        slack = find(bus[:, BUS_TYPE] == REF)
        slack = slack[0]

    ## set the slack bus to be used to compute initial PTDF
    if isscalar(slack):
        slack_bus = slack
    else:
        slack_bus = 0      ## use bus 1 for temp slack bus

    nb = bus.shape[0]
    nbr = branch.shape[0]
    noref = arange(1, nb)      ## use bus 1 for voltage angle reference
    noslack = find(arange(nb) != slack_bus)

    ## check that bus numbers are equal to indices to bus (one set of bus numbers)
    if any(bus[:, BUS_I] != arange(nb)):
        stderr.write('makePTDF: buses must be numbered consecutively')

    ## compute PTDF for single slack_bus
    Bbus, Bf, _, _ = makeBdc(baseMVA, bus, branch)
    Bbus, Bf = Bbus.todense(), Bf.todense()
    H = zeros((nbr, nb))
    H[:, noslack] = solve( Bbus[ix_(noslack, noref)].T, Bf[:, noref].T ).T
    #             = Bf[:, noref] * inv(Bbus[ix_(noslack, noref)])

    ## distribute slack, if requested
    if not isscalar(slack):
        if len(slack.shape) == 1:  ## slack is a vector of weights
            slack = slack / sum(slack)   ## normalize weights

            ## conceptually, we want to do ...
            ##    H = H * (eye(nb, nb) - slack * ones((1, nb)))
            ## ... we just do it more efficiently
            v = dot(H, slack)
            for k in range(nb):
                H[:, k] = H[:, k] - v
        else:
            H = dot(H, slack)

    return H

    def _sys_event(self, C):
        sysdef = self.sysdef
        row, col = C.shape

        if self.systype.lower() == "series":
            self.cutset = 1 - np.prod(np.ones((row, col)) - C, axis=1, dtype=int)
            ncutsets = []
        elif self.systype.lower() == "parallel":
            self.cutset = np.prod(C, axis=1, dtype=int)
            ncutsets = []
        else:
            if self.sysdef[1].lower() == "link":
                self.sysdef = -self.sysdef
            self.sysdef[0] = np.hstack((0, self.sysdef[0], 0))
            sysnonzero = np.find(self.sysdef[0] != 0)[0]
            syszero = np.find(self.sysdef[0] == 0)[0]
            int1 = syszero - np.hstack((0, syszero[:-1]))
            sizeCutSets = int1[int1 > 1] - 1
            ncutsets = sizeCutSets.shape[0]
            for i in xrange(ncutsets):
                cCutSet = np.ones(row, 1)
                for j in xrange(sizeCutSets[i]):
                    comp = self.sysdef[0][sysnonzero[np.sum(sizeCutSets[:i]) + j]]
                    if comp < 0:
                        cCutSet = cCutSet * (np.ones((row, 1)) - C.T[:, abs(comp)])
                    else:
                        cCutSet = cCutSet * C.T[:, comp]
                cCutSets[:, i] = cCutSet
            self.cutset = np.ones(row, 1) - np.prod(np.ones(row, ncutsets) - cCutSets, axis=1)

def rhoa( snd, cal=260e-9, corrramp=True ):
    ''' compute apparent resistivity from sounding '''
    Tx = np.prod( [float(a) for a in snd['LOOP_SIZE'].split()] )
    if 'COIL_SIZE' in snd:
        Rx = snd['COIL_SIZE']
    else:
        Rx = Tx
    
    v = snd['VOLTAGE']
    istart, istop = 0, len(v) # default: take all
    mav = np.find( v == max(v) )
    if len(mav) > 1: #several equal big ones: start after
        istart = max(mav)+1

    if min(v) < 0.0: # negative values: stop at first
        istop = min( np.find( v < 0.0 ) )
    
    v = v[istart:istop]
    dv = snd['ST_DEV'][istart:istop] #/ snd['CURRENT']
    t = snd['TIME'][istart:istop]
    if corrramp and 'RAMP_TIME' in snd: 
        t = t - snd['RAMP_TIME']
    
    if Rx == 1: # apparently B-field
        rhoa = rhoafromB( v * cal, t, Tx )
    else:
        rhoa = rhoafromU( v, t, Tx, Rx )

    rhoaerr = dv / v * (2./3.)
    return rhoa, t, rhoaerr

def bustypes(bus, gen):
    """Builds index lists of each type of bus (C{REF}, C{PV}, C{PQ}).

    Generators with "out-of-service" status are treated as L{PQ} buses with
    zero generation (regardless of C{Pg}/C{Qg} values in gen). Expects C{bus}
    and C{gen} have been converted to use internal consecutive bus numbering.

    @param bus: bus data
    @param gen: generator data
    @return: index lists of each bus type

    @author: Ray Zimmerman (PSERC Cornell)
    """
    # get generator status
    nb = bus.shape[0]
    ng = gen.shape[0]
    # gen connection matrix, element i, j is 1 if, generator j at bus i is ON
    Cg = sparse((gen[:, GEN_STATUS] > 0,
                 (gen[:, GEN_BUS], range(ng))), (nb, ng))
    # number of generators at each bus that are ON
    bus_gen_status = (Cg * ones(ng, int)).astype(bool)

    # form index lists for slack, PV, and PQ buses
    ref = find((bus[:, BUS_TYPE] == REF) & bus_gen_status) # ref bus index
    pv  = find((bus[:, BUS_TYPE] == PV)  & bus_gen_status) # PV bus indices
    pq  = find((bus[:, BUS_TYPE] == PQ) | ~bus_gen_status) # PQ bus indices

    # pick a new reference bus if for some reason there is none (may have been
    # shut down)
    if len(ref) == 0:
        ref = pv[0]      # use the first PV bus
        pv = pv[1:]      # take it off PV list

    return ref, pv, pq

    def __init__(self, paramfile):
        for row in paramfile:

            if row[0]=='BIAS_COMBINE': self.BIAS_COMBINE = row[1]
            if row[0]=='COMP_COMBINE': self.COMP_COMBINE = row[1]
            if row[0]=='FLAT_COMBINE': self.FLAT_COMBINE = row[1]
            if row[0]=='FIBER_TRACE': self.FIBER_TRACE = row[1]
            if row[0]=='BIAS_SUBTRACT': self.BIAS_SUBTRACT = row[1]
            if row[0]=='LAMBDA_SOLVE': self.LAMBDA_SOLVE = row[1]
            if row[0]=='COSMIC_RAYS': self.COSMIC_RAYS = row[1]
            if row[0]=='SKY_SUBTRACT': self.SKY_SUBTRACT = row[1]
            if row[0]=='WRITE_SPECTRA': self.WRITE_SPECTRA = row[1]
            
            if row[0]=='FIBERS_TOTAL': self.FIBERS_TOTAL = int(row[1])
            if row[0]=='FIBERS_EXCLUDE':
                fibs = np.array( [ int(i) for i in row[1].split(',') ] )
                fibs.sort()
                fibs = fibs[ np.find( 0<fibs ) ]
                fibs = fibs[ np.find( fibs<self.FIBERS_TOTAL+1 ) ]
                self.FIBERS_EXCLUDE = fibs
            if row[0]=='FIBER_WIDTH': self.FIBER_WIDTH = int(row[1])
            if row[0]=='FIBER_SEP': self.FIBER_SEP = int(row[1])
            if row[0]=='LO_BUFFER': self.LO_BUFFER = int(row[1])
            if row[0]=='HI_BUFFER': self.HI_BUFFER = int(row[1])
            if row[0]=='LINE_LIST': self.LINE_LIST = row[1]
            if row[0]=='COMP_HEADER': self.COMP_HEADER = row[1]
            if row[0]=='COMP_NAME': self.COMP_NAME = row[1]
            if row[0]=='POLYFIT_ORDER': self.POLYFIT_ORDER = int(row[1])
            if row[0]=='CR_SCAN_DX': self.CR_SCAN_DX = int(row[1])

def bustypes(bus, gen, Sbus):
    """
    Builds index lists of each type of bus (C{REF}, C{PV}, C{PQ}).

    Generators with "out-of-service" status are treated as L{PQ} buses with
    zero generation (regardless of C{Pg}/C{Qg} values in gen). Expects C{bus}
    and C{gen} have been converted to use internal consecutive bus numbering.

    @param bus: bus data
    @param gen: generator data
    @return: index lists of each bus type

    @author: Ray Zimmerman (PSERC Cornell)
    """
    # flag to indicate that it is impossible to solve the grid
    the_grid_is_disabled = False

    # get generator status
    nb = bus.shape[0]
    ng = gen.shape[0]
    # gen connection matrix, element i, j is 1 if, generator j at bus i is ON
    Cg = sparse((gen[:, GEN_STATUS] > 0,
                 (gen[:, GEN_BUS], range(ng))), (nb, ng))
    # number of generators at each bus that are ON
    bus_gen_status = (Cg * ones(ng, int)).astype(bool)

    # form index lists for slack, PV, and PQ buses
    ref = find((bus[:, BUS_TYPE] == REF) & bus_gen_status)  # ref bus index
    pv = find((bus[:, BUS_TYPE] == PV) & bus_gen_status)  # PV bus indices
    pq = find((bus[:, BUS_TYPE] == PQ) | ~bus_gen_status)  # PQ bus indices

    # pick a new reference bus if for some reason there is none (may have been
    # shut down)
    if len(ref) == 0:
        if len(pv) > 0:
            ref = [pv[0]]      # use the first PV bus
            pv = pv[1:]      # take it off PV list
        else:   # look for positive power injections to take the largest as the slack
            positive_power_injections = Sbus.real[where(Sbus.real > 0)[0]]
            if len(positive_power_injections) > 0:
                idx = where(Sbus.real == max(positive_power_injections))[0]
                if len(idx) == 1:
                    ref = idx
                    i = where(pq == idx[0])[0][0]
                    pq = delete(pq, i)
                else:
                    warn('It was not possible to find a slack bus')
                    the_grid_is_disabled = True
            else:
                warn('It was not possible to find a slack bus')
                the_grid_is_disabled = True

    # create the types array
    types = zeros(nb)
    types[ref] = 3
    types[pv] = 2
    types[pq] = 1

    return ref, pv, pq, types, the_grid_is_disabled

def userfcn_reserves_ext2int(ppc, *args):
    """This is the 'ext2int' stage userfcn callback that prepares the input
    data for the formulation stage. It expects to find a 'reserves' field
    in ppc as described above. The optional args are not currently used.
    """
    ## initialize some things
    r = ppc['reserves']
    o = ppc['order']
    ng0 = o['ext']['gen'].shape[0]  ## number of original gens (+ disp loads)
    nrz = r['req'].shape[0]         ## number of reserve zones
    if nrz > 1:
        ppc['reserves']['rgens'] = any(r['zones'], 0)  ## mask of gens available to provide reserves
    else:
        ppc['reserves']['rgens'] = r['zones']

    igr = find(ppc['reserves']['rgens']) ## indices of gens available to provide reserves
    ngr = len(igr)              ## number of gens available to provide reserves

    ## check data for consistent dimensions
    if r['zones'].shape[0] != nrz:
        stderr.write('userfcn_reserves_ext2int: the number of rows in ppc[\'reserves\'][\'req\'] (%d) and ppc[\'reserves\'][\'zones\'] (%d) must match\n' % (nrz, r['zones'].shape[0]))

    if (r['cost'].shape[0] != ng0) & (r['cost'].shape[0] != ngr):
        stderr.write('userfcn_reserves_ext2int: the number of rows in ppc[\'reserves\'][\'cost\'] (%d) must equal the total number of generators (%d) or the number of generators able to provide reserves (%d)\n' % (r['cost'].shape[0], ng0, ngr))

    if 'qty' in r:
        if r['qty'].shape[0] != r['cost'].shape[0]:
            stderr.write('userfcn_reserves_ext2int: ppc[\'reserves\'][\'cost\'] (%d x 1) and ppc[\'reserves\'][\'qty\'] (%d x 1) must be the same dimension\n' % (r['cost'].shape[0], r['qty'].shape[0]))


    ## convert both cost and qty from ngr x 1 to full ng x 1 vectors if necessary
    if r['cost'].shape[0] < ng0:
        if 'original' not in ppc['reserves']:
            ppc['reserves']['original'] = {}
        ppc['reserves']['original']['cost'] = r['cost'].copy()    ## save original
        cost = zeros(ng0)
        cost[igr] = r['cost']
        ppc['reserves']['cost'] = cost
        if 'qty' in r:
            ppc['reserves']['original']['qty'] = r['qty'].copy()  ## save original
            qty = zeros(ng0)
            qty[igr] = r['qty']
            ppc['reserves']['qty'] = qty

    ##-----  convert stuff to internal indexing  -----
    ## convert all reserve parameters (zones, costs, qty, rgens)
    if 'qty' in r:
        ppc = e2i_field(ppc, ['reserves', 'qty'], 'gen')

    ppc = e2i_field(ppc, ['reserves', 'cost'], 'gen')
    ppc = e2i_field(ppc, ['reserves', 'zones'], 'gen', 1)
    ppc = e2i_field(ppc, ['reserves', 'rgens'], 'gen', 1)

    ## save indices of gens available to provide reserves
    ppc['order']['ext']['reserves']['igr'] = igr               ## external indexing
    ppc['reserves']['igr'] = find(ppc['reserves']['rgens'])    ## internal indexing

    return ppc

def modcost(gencost, alpha, modtype='SCALE_F'):
    """Modifies generator costs by shifting or scaling (F or X).

    For each generator cost F(X) (for real or reactive power) in
    C{gencost}, this function modifies the cost by scaling or shifting
    the function by C{alpha}, depending on the value of C{modtype}, and
    and returns the modified C{gencost}. Rows of C{gencost} can be a mix
    of polynomial or piecewise linear costs.

    C{modtype} takes one of the 4 possible values (let F_alpha(X) denote the
    the modified function)::
        SCALE_F (default) : F_alpha(X)         == F(X) * ALPHA
        SCALE_X           : F_alpha(X * ALPHA) == F(X)
        SHIFT_F           : F_alpha(X)         == F(X) + ALPHA
        SHIFT_X           : F_alpha(X + ALPHA) == F(X)

    @author: Ray Zimmerman (PSERC Cornell)
    @author: Richard Lincoln
    """
    gencost = gencost.copy()

    ng, m = gencost.shape
    if ng != 0:
        ipwl = find(gencost[:, MODEL] == PW_LINEAR)
        ipol = find(gencost[:, MODEL] == POLYNOMIAL)
        c = gencost[ipol, COST:m]

        if modtype == 'SCALE_F':
            gencost[ipol, COST:m]       = alpha * c
            gencost[ipwl, COST+1:m:2]   = alpha * gencost[ipwl, COST + 1:m:2]
        elif modtype == 'SCALE_X':
            for k in range(len(ipol)):
                n = gencost[ipol[k], NCOST].astype(int)
                for i in range(n):
                    gencost[ipol[k], COST + i] = c[k, i] / alpha**(n - i - 1)
            gencost[ipwl, COST:m - 1:2]   = alpha * gencost[ipwl, COST:m - 1:2]
        elif modtype == 'SHIFT_F':
            for k in range(len(ipol)):
                n = gencost[ipol[k], NCOST].astype(int)
                gencost[ipol[k], COST + n - 1] = alpha + c[k, n - 1]
            gencost[ipwl, COST+1:m:2]   = alpha + gencost[ipwl, COST + 1:m:2]
        elif modtype == 'SHIFT_X':
            for k in range(len(ipol)):
                n = gencost[ipol[k], NCOST].astype(int)
                gencost[ipol[k], COST:COST + n] = \
                        polyshift(c[k, :n].T, alpha).T
            gencost[ipwl, COST:m - 1:2]   = alpha + gencost[ipwl, COST:m - 1:2]
        else:
            sys.stderr.write('modcost: "%s" is not a valid modtype\n' % modtype)

    return gencost

def makeSbus(baseMVA, bus, gen):
    """Builds the vector of complex bus power injections.

    Returns the vector of complex bus power injections, that is, generation
    minus load. Power is expressed in per unit.

    @see: L{makeYbus}

    @author: Ray Zimmerman (PSERC Cornell)
    @author: Richard Lincoln
    """
    ## generator info
    on = find(gen[:, GEN_STATUS] > 0)      ## which generators are on?
    gbus = gen[on, GEN_BUS]                   ## what buses are they at?

    ## form net complex bus power injection vector
    nb = bus.shape[0]
    ngon = on.shape[0]
    ## connection matrix, element i, j is 1 if gen on(j) at bus i is ON
    Cg = sparse((ones(ngon), (gbus, range(ngon))), (nb, ngon))

    ## power injected by gens plus power injected by loads converted to p.u.
    Sbus = ( Cg * (gen[on, PG] + 1j * gen[on, QG]) -
             (bus[:, PD] + 1j * bus[:, QD]) ) / baseMVA

    return Sbus

def userfcn_iflims_ext2int(ppc, *args):
    """This is the 'ext2int' stage userfcn callback that prepares the input
    data for the formulation stage. It expects to find an 'if' field in
    ppc as described above. The optional args are not currently used.
    """
    ## initialize some things
    ifmap = ppc['if']['map']
    o = ppc['order']
    nl0 = o['ext']['branch'].shape[0]    ## original number of branches
    nl = ppc['branch'].shape[0]          ## number of on-line branches

    ## save if.map for external indexing
    ppc['order']['ext']['ifmap'] = ifmap

    ##-----  convert stuff to internal indexing  -----
    e2i = zeros(nl0)
    e2i[o['branch']['status']['on']] = arange(nl)  ## ext->int branch index mapping
    d = sign(ifmap[:, 1])
    br = abs(ifmap[:, 1]).astype(int)
    ifmap[:, 1] = d * e2i[br]

    ifmap = delete(ifmap, find(ifmap[:, 1] == 0), 0)  ## delete branches that are out

    ppc['if']['map'] = ifmap

    return ppc

def userfcn_reserves_formulation(om, *args):
    """This is the 'formulation' stage userfcn callback that defines the
    user costs and constraints for fixed reserves. It expects to find
    a 'reserves' field in the ppc stored in om, as described above.
    By the time it is passed to this callback, ppc['reserves'] should
    have two additional fields:
        - C{igr}     C{1 x ngr}, indices of generators available for reserves
        - C{rgens}   C{1 x ng}, 1 if gen avaiable for reserves, 0 otherwise
    It is also assumed that if cost or qty were C{ngr x 1}, they have been
    expanded to C{ng x 1} and that everything has been converted to
    internal indexing, i.e. all gens are on-line (by the 'ext2int'
    callback). The optional args are not currently used.
    """
    ## initialize some things
    ppc = om.get_ppc()
    r = ppc['reserves']
    igr = r['igr']                ## indices of gens available to provide reserves
    ngr = len(igr)                ## number of gens available to provide reserves
    ng  = ppc['gen'].shape[0]     ## number of on-line gens (+ disp loads)

    ## variable bounds
    Rmin = zeros(ngr)               ## bound below by 0
    Rmax = Inf * ones(ngr)          ## bound above by ...
    k = find(ppc['gen'][igr, RAMP_10])
    Rmax[k] = ppc['gen'][igr[k], RAMP_10] ## ... ramp rate and ...
    if 'qty' in r:
        k = find(r['qty'][igr] < Rmax)
        Rmax[k] = r['qty'][igr[k]]        ## ... stated max reserve qty
    Rmax = Rmax / ppc['baseMVA']

    ## constraints
    I = speye(ngr, ngr, format='csr')                     ## identity matrix
    Ar = hstack([sparse((ones(ngr), (arange(ngr), igr)), (ngr, ng)), I], 'csr')
    ur = ppc['gen'][igr, PMAX] / ppc['baseMVA']
    lreq = r['req'] / ppc['baseMVA']

    ## cost
    Cw = r['cost'][igr] * ppc['baseMVA']     ## per unit cost coefficients

    ## add them to the model
    om.add_vars('R', ngr, [], Rmin, Rmax)
    om.add_constraints('Pg_plus_R', Ar, [], ur, ['Pg', 'R'])
    om.add_constraints('Rreq', sparse( r['zones'][:, igr] ), lreq, [], ['R'])
    om.add_costs('Rcost', {'N': I, 'Cw': Cw}, ['R'])

    return om

def makeAang(baseMVA, branch, nb, ppopt):
    """Construct constraints for branch angle difference limits.

    Constructs the parameters for the following linear constraint limiting
    the voltage angle differences across branches, where C{Va} is the vector
    of bus voltage angles. C{nb} is the number of buses::

        lang <= Aang * Va <= uang

    C{iang} is the vector of indices of branches with angle difference limits.

    @author: Ray Zimmerman (PSERC Cornell)
    @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad
    Autonoma de Manizales)
    @author: Richard Lincoln
    """
    ## options
    ignore_ang_lim = ppopt['OPF_IGNORE_ANG_LIM']

    if ignore_ang_lim:
        Aang  = zeros((0, nb))
        lang  = array([])
        uang  = array([])
        iang  = array([])
    else:
        iang = find(((branch[:, ANGMIN] != 0) & (branch[:, ANGMIN] > -360)) |
                    ((branch[:, ANGMAX] != 0) & (branch[:, ANGMAX] <  360)))
        iangl = find(branch[iang, ANGMIN])
        iangh = find(branch[iang, ANGMAX])
        nang = len(iang)

        if nang > 0:
            ii = r_[arange(nang), arange(nang)]
            jj = r_[branch[iang, F_BUS], branch[iang, T_BUS]]
            Aang = sparse((r_[ones(nang), -ones(nang)],
                           (ii, jj)), (nang, nb))
            uang = Inf * ones(nang)
            lang = -uang
            lang[iangl] = branch[iang[iangl], ANGMIN] * pi / 180
            uang[iangh] = branch[iang[iangh], ANGMAX] * pi / 180
        else:
            Aang  = zeros((0, nb))
            lang  = array([])
            uang  = array([])

    return Aang, lang, uang, iang

def userfcn_dcline_int2ext(results, args):
    """This is the 'int2ext' stage userfcn callback that converts everything
    back to external indexing and packages up the results. It expects to
    find a 'dcline' field in the results struct as described for ppc
    above. It also expects that the last 2*ndc entries in the gen and
    gencost matrices correspond to the in-service DC lines (where ndc is
    the number of rows in MPC.dcline. These extra rows are removed from
    gen and gencost and the flow is taken from the PG of these gens and
    placed in the flow column of the appropiate dcline row. The
    optional args are not currently used.
    """
    c = idx_dcline.c

    ## initialize some things
    o = results['order']
    k = find(o['ext']['dcline'][:, c['BR_STATUS']])
    ndc = len(k)                    ## number of in-service DC lines
    ng  = results['gen'].shape[0] - 2*ndc; ## number of original gens/disp loads

    ## extract dummy gens
    fg = results['gen'][ng:ng + ndc, :]
    tg = results['gen'][ng + ndc:ng + 2 * ndc, :]

    ## remove dummy gens
    #results['gen']     = results['gen'][:ng + 1, :]
    #results['gencost'] = results['gencost'][:ng + 1, :]
    results['gen']     = results['gen'][:ng, :]
    results['gencost'] = results['gencost'][:ng, :]

    ## get the solved flows
    results['dcline'][:, c['PF']] = -fg[:, PG]
    results['dcline'][:, c['PT']] =  tg[:, PG]
    results['dcline'][:, c['QF']] =  fg[:, QG]
    results['dcline'][:, c['QT']] =  tg[:, QG]
    results['dcline'][:, c['VF']] =  fg[:, VG]
    results['dcline'][:, c['VT']] =  tg[:, VG]
    if fg.shape[1] >= MU_QMIN:
        results['dcline'] = c_[results['dcline'], zeros((ndc, 6))]
        results['dcline'][:, c['MU_PMIN'] ] = fg[:, MU_PMAX] + tg[:, MU_PMIN]
        results['dcline'][:, c['MU_PMAX'] ] = fg[:, MU_PMIN] + tg[:, MU_PMAX]
        results['dcline'][:, c['MU_QMINF']] = fg[:, MU_QMIN]
        results['dcline'][:, c['MU_QMAXF']] = fg[:, MU_QMAX]
        results['dcline'][:, c['MU_QMINT']] = tg[:, MU_QMIN]
        results['dcline'][:, c['MU_QMAXT']] = tg[:, MU_QMAX]

    results['order']['int'] = {}
    ##-----  convert stuff back to external indexing  -----
    results['order']['int']['dcline'] = results['dcline']  ## save internal version
    ## copy results to external version
    o['ext']['dcline'][k, c['PF']:c['VT'] + 1] = results['dcline'][:, c['PF']:c['VT'] + 1]
    if results['dcline'].shape[1] == c['MU_QMAXT'] + 1:
        o['ext']['dcline'] = c_[o['ext']['dcline'], zeros((ndc, 6))]
        o['ext']['dcline'][k, c['MU_PMIN']:c['MU_QMAXT'] + 1] = \
                results['dcline'][:, c['MU_PMIN']:c['MU_QMAXT'] + 1]

    results['dcline'] = o['ext']['dcline']            ## use external version

    return results

def totcost(gencost, Pg):
    """Computes total cost for generators at given output level.

    Computes total cost for generators given a matrix in gencost format and
    a column vector or matrix of generation levels. The return value has the
    same dimensions as PG. Each row of C{gencost} is used to evaluate the
    cost at the points specified in the corresponding row of C{Pg}.

    @author: Ray Zimmerman (PSERC Cornell)
    @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad
    Autonoma de Manizales)
    @author: Richard Lincoln
    """
    ng, m = gencost.shape
    totalcost = zeros(ng)

    if len(gencost) > 0:
        ipwl = find(gencost[:, MODEL] == PW_LINEAR)
        ipol = find(gencost[:, MODEL] == POLYNOMIAL)
        if len(ipwl) > 0:
            p = gencost[:, COST:(m-1):2]
            c = gencost[:, (COST+1):m:2]

            for i in ipwl:
                ncost = gencost[i, NCOST]
                for k in arange(ncost - 1):
                    p1, p2 = p[i, k], p[i, k+1]
                    c1, c2 = c[i, k], c[i, k+1]
                    m = (c2 - c1) / (p2 - p1)
                    b = c1 - m * p1
                    Pgen = Pg[i]
                    if Pgen < p2:
                        totalcost[i] = m * Pgen + b
                        break
                    totalcost[i] = m * Pgen + b

        if len(ipol) > 0:
            totalcost[ipol] = polycost(gencost[ipol, :], Pg[ipol])

    return totalcost

def polycost(gencost, Pg, der=0):
    """Evaluates polynomial generator cost & derivatives.

    C{f = polycost(gencost, Pg)} returns the vector of costs evaluated at C{Pg}

    C{df = polycost(gencost, Pg, 1)} returns the vector of first derivatives
    of costs evaluated at C{Pg}

    C{d2f = polycost(gencost, Pg, 2)} returns the vector of second derivatives
    of costs evaluated at C{Pg}

    C{gencost} must contain only polynomial costs
    C{Pg} is in MW, not p.u. (works for C{Qg} too)

    @author: Ray Zimmerman (PSERC Cornell)
    @author: Richard Lincoln
    """
    if any(gencost[:, MODEL] == PW_LINEAR):
        sys.stderr.write('polycost: all costs must be polynomial\n')

    ng = len(Pg)
    maxN = max( gencost[:, NCOST].astype(int) )
    minN = min( gencost[:, NCOST].astype(int) )

    ## form coefficient matrix where 1st column is constant term, 2nd linear, etc.
    c = zeros((ng, maxN))
    for n in arange(minN, maxN + 1):
        k = find(gencost[:, NCOST] == n)   ## cost with n coefficients
        c[k, :n] = gencost[k, (COST + n - 1):COST - 1:-1]

    ## do derivatives
    for d in range(1, der + 1):
        if c.shape[1] >= 2:
            c = c[:, 1:maxN - d + 1]
        else:
            c = zeros((ng, 1))
            break

        for k in range(2, maxN - d + 1):
            c[:, k-1] = c[:, k-1] * k

    ## evaluate polynomial
    if len(c) == 0:
        f = zeros(Pg.shape)
    else:
        f = c[:, :1].flatten()  ## constant term
        for k in range(1, c.shape[1]):
            f = f + c[:, k] * Pg**k

    return f