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

def despike(datax,datay,threshold=3.0,interpolate=False):
    """
    Remove spikes larger than threshold*(standard deviation of diffs in datay).
    interpolate=False (default) just removes the points;
    interpolate=True will just make them equal to the point to the left 
    (if you need to keep array length constant, say).
    """
    d = np.diff(datay)
    spikes = np.find(abs(d)>threshold*std(d)) 
    if (len(spikes)>0.1*len(d)): 
        print 'Did not remove spikes because it wanted to remove too many.'
        return datax,datay # if we're trying to remove a lot of points, just send it back...
    spikes=np.delete(spikes,find(diff(spikes)==1)+1) # if spike is one point, don't delete point after as well.
    if interpolate==False:
        datax = np.delete(datax,spikes+1)
        datay = np.delete(datay,spikes+1)
    else:
        # don't actually 'interpolate', just set it equal to the previous point
        # (actual interpolation could get messy in the case of a two-point spike, for example)
        for i in spikes+1: datay[i] = datay[i-3]
    if (np.any(np.abs(diff(datay))>threshold*np.std(d))): datax,datay = despike(datax,datay,threshold,interpolate)
    return datax, datay

#.........这里部分代码省略.........

    t = ''.join([t0, 'AC PF (with DC lines) : '])
    ppc1 = {'baseMVA': r['baseMVA'],
            'bus': r['bus'][:, :VMIN + 1].copy(),
            'gen': r['gen'][:, :APF + 1].copy(),
            'branch': r['branch'][:, :ANGMAX + 1].copy(),
            'gencost': r['gencost'].copy(),
            'dcline': r['dcline'][:, :c.LOSS1 + 1].copy()}
    ppc1 = toggle_dcline(ppc1, 'on')
    ppc1['bus'][:, VM] = 1
    ppc1['bus'][:, VA] = 0
    rp = runpf(ppc1, ppopt)
    success = rp['success']
    t_ok(success, [t, 'success'])
    t_is(   rp['bus'][:,ib_voltage],    r['bus'][:,ib_voltage], 3, [t, 'bus voltage'])
    #t_is(   rp['gen'][:,ig_disp   ],    r['gen'][:,ig_disp   ], 3, [t, 'gen dispatch'])
    t_is(   rp['gen'][:2,ig_disp ],    r['gen'][:2,ig_disp ], 3, [t, 'gen dispatch'])
    t_is(   rp['gen'][2,PG        ],    r['gen'][2,PG        ], 3, [t, 'gen dispatch'])
    t_is(   rp['gen'][2,QG]+rp['dcline'][0,c.QF], r['gen'][2,QG]+r['dcline'][0,c.QF], 3, [t, 'gen dispatch'])
    t_is(rp['branch'][:,ibr_flow  ], r['branch'][:,ibr_flow  ], 3, [t, 'branch flow'])

    ## add appropriate P and Q injections and check angles and generation when running PF
    t = ''.join([t0, 'AC PF (with equivalent injections) : '])
    ppc1 = {'baseMVA': r['baseMVA'],
            'bus': r['bus'][:, :VMIN + 1].copy(),
            'gen': r['gen'][:, :APF + 1].copy(),
            'branch': r['branch'][:, :ANGMAX + 1].copy(),
            'gencost': r['gencost'].copy(),
            'dcline': r['dcline'][:, :c.LOSS1 + 1].copy()}
    ppc1['bus'][:, VM] = 1
    ppc1['bus'][:, VA] = 0
    for k in range(ndc):
        if ppc1['dcline'][k, c.BR_STATUS]:
            ff = find(ppc1['bus'][:, BUS_I] == ppc1['dcline'][k, c.F_BUS])
            tt = find(ppc1['bus'][:, BUS_I] == ppc1['dcline'][k, c.T_BUS])
            ppc1['bus'][ff, PD] = ppc1['bus'][ff, PD] + r['dcline'][k, c.PF]
            ppc1['bus'][ff, QD] = ppc1['bus'][ff, QD] - r['dcline'][k, c.QF]
            ppc1['bus'][tt, PD] = ppc1['bus'][tt, PD] - r['dcline'][k, c.PT]
            ppc1['bus'][tt, QD] = ppc1['bus'][tt, QD] - r['dcline'][k, c.QT]
            ppc1['bus'][ff, VM] = r['dcline'][k, c.VF]
            ppc1['bus'][tt, VM] = r['dcline'][k, c.VT]
            ppc1['bus'][ff, BUS_TYPE] = PV
            ppc1['bus'][tt, BUS_TYPE] = PV

    rp = runpf(ppc1, ppopt)
    success = rp['success']
    t_ok(success, [t, 'success'])
    t_is(   rp['bus'][:,ib_voltage],    r['bus'][:,ib_voltage],  3, [t, 'bus voltage'])
    t_is(   rp['gen'][:,ig_disp   ],    r['gen'][:,ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(rp['branch'][:,ibr_flow  ], r['branch'][:,ibr_flow  ],  3, [t, 'branch flow'])

    ## test DC OPF
    t = ''.join([t0, 'DC OPF (with DC lines) : '])
    ppc = ppc0.copy()
    ppc['gen'][0, PMIN] = 10
    ppc['branch'][4, RATE_A] = 100
    ppc = toggle_dcline(ppc, 'on')
    r = rundcopf(ppc, ppopt)
    success = r['success']
    t_ok(success, [t, 'success'])
    expected = array([
        [10, 8.9, 0, 0, 1.01, 1],
        [2,  2,   0, 0, 1,    1],
        [0,  0,   0, 0, 1,    1],
        [10, 9.5, 0, 0, 1, 0.98]
    ])

#.........这里部分代码省略.........

        if len(l) == 0:
            l = -Inf * ones(nA)     ## ... unbounded below.

    if len(x0) == 0:
        x0 = zeros(nx)

    ## default options
    if 'verbose' in opt:
        verbose = opt['verbose']
    else:
        verbose = 0

    if 'max_it' in opt:
        max_it = opt['max_it']
    else:
        max_it = 0

    ## make sure args are sparse/full as expected by IPOPT
    if len(H) > 0:
        if not issparse(H):
            H = sparse(H)

    if not issparse(A):
        A = sparse(A)

    ##-----  run optimization  -----
    ## set options dict for IPOPT
    options = {}
    if 'ipopt_opt' in opt:
        options['ipopt'] = ipopt_options(opt['ipopt_opt'])
    else:
        options['ipopt'] = ipopt_options()

    options['ipopt']['jac_c_constant']    = 'yes'
    options['ipopt']['jac_d_constant']    = 'yes'
    options['ipopt']['hessian_constant']  = 'yes'
    options['ipopt']['least_square_init_primal']  = 'yes'
    options['ipopt']['least_square_init_duals']   = 'yes'
    # options['ipopt']['mehrotra_algorithm']        = 'yes'     ## default 'no'
    if verbose:
        options['ipopt']['print_level'] = min(12, verbose * 2 + 1)
    else:
        options['ipopt']['print_level = 0']

    if max_it:
        options['ipopt']['max_iter'] = max_it

    ## define variable and constraint bounds, if given
    if nA:
        options['cu'] = u
        options['cl'] = l

    if len(xmin) > 0:
        options['lb'] = xmin

    if len(xmax) > 0:
        options['ub'] = xmax

    ## assign function handles
    funcs = {}
    funcs['objective']         = lambda x: 0.5 * x.T * H * x + c.T * x
    funcs['gradient']          = lambda x: H * x + c
    funcs['constraints']       = lambda x: A * x
    funcs['jacobian']          = lambda x: A
    funcs['jacobianstructure'] = lambda : A
    funcs['hessian']           = lambda x, sigma, lmbda: tril(H)
    funcs['hessianstructure']  = lambda : tril(H)

    ## run the optimization
    x, info = pyipopt(x0, funcs, options)

    if info['status'] == 0 | info['status'] == 1:
        eflag = 1
    else:
        eflag = 0

    output = {}
    if 'iter' in info:
        output['iterations'] = info['iter']

    output['info']       = info['status']
    f = funcs['objective'](x)

    ## repackage lmbdas
    kl = find(info['lmbda'] < 0)                     ## lower bound binding
    ku = find(info['lmbda'] > 0)                     ## upper bound binding
    mu_l = zeros(nA)
    mu_l[kl] = -info['lmbda'][kl]
    mu_u = zeros(nA)
    mu_u[ku] = info['lmbda'][ku]

    lmbda = {
        'mu_l': mu_l,
        'mu_u': mu_u,
        'lower': info['zl'],
        'upper': info['zu']
    }

    return x, f, eflag, output, lmbda

def pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq):
    """Updates bus, gen, branch data structures to match power flow soln.

    @author: Ray Zimmerman (PSERC Cornell)
    @author: Richard Lincoln
    """
    ## initialize return values
    bus     = bus0
    gen     = gen0
    branch  = branch0

    ##----- update bus voltages -----
    bus[:, VM] = abs(V)
    bus[:, VA] = angle(V) * 180 / pi

    ##----- update Qg for all gens and Pg for slack bus(es) -----
    ## generator info
    on = find(gen[:, GEN_STATUS] > 0) ## which generators are on?
    gbus = gen[on, GEN_BUS].astype(int)  ## what buses are they at?

    ## compute total injected bus powers
    Sbus = V[gbus] * conj(Ybus[gbus, :] * V)

    ## update Qg for all generators
    gen[:, QG] = zeros(gen.shape[0])              ## zero out all Qg
    gen[on, QG] = Sbus.imag * baseMVA + bus[gbus, QD]    ## inj Q + local Qd
    ## ... at this point any buses with more than one generator will have
    ## the total Q dispatch for the bus assigned to each generator. This
    ## must be split between them. We do it first equally, then in proportion
    ## to the reactive range of the generator.

    if len(on) > 1:
        ## build connection matrix, element i, j is 1 if gen on(i) at bus j is ON
        nb = bus.shape[0]
        ngon = on.shape[0]
        Cg = csr_matrix((ones(ngon), (range(ngon), gbus)), (ngon, nb))

        ## divide Qg by number of generators at the bus to distribute equally
        ngg = Cg * Cg.sum(0).T    ## ngon x 1, number of gens at this gen's bus
        ngg = asarray(ngg).flatten()  # 1D array
        gen[on, QG] = gen[on, QG] / ngg

        ## divide proportionally
        Cmin = csr_matrix((gen[on, QMIN], (range(ngon), gbus)), (ngon, nb))
        Cmax = csr_matrix((gen[on, QMAX], (range(ngon), gbus)), (ngon, nb))
        Qg_tot = Cg.T * gen[on, QG]## nb x 1 vector of total Qg at each bus
        Qg_min = Cmin.sum(0).T       ## nb x 1 vector of min total Qg at each bus
        Qg_max = Cmax.sum(0).T       ## nb x 1 vector of max total Qg at each bus
        Qg_min = asarray(Qg_min).flatten()  # 1D array
        Qg_max = asarray(Qg_max).flatten()  # 1D array
        ## gens at buses with Qg range = 0
        ig = find(Cg * Qg_min == Cg * Qg_max)
        Qg_save = gen[on[ig], QG]
        gen[on, QG] = gen[on, QMIN] + \
            (Cg * ((Qg_tot - Qg_min) / (Qg_max - Qg_min + EPS))) * \
                (gen[on, QMAX] - gen[on, QMIN])    ##    ^ avoid div by 0
        gen[on[ig], QG] = Qg_save  ## (terms are mult by 0 anyway)

    ## update Pg for slack bus(es)
    ## inj P + local Pd
    for k in range(len(ref)):
        refgen = find(gbus == ref[k])  ## which is(are) the reference gen(s)?
        gen[on[refgen[0]], PG] = \
                Sbus[refgen[0]].real * baseMVA + bus[ref[k], PD]
        if len(refgen) > 1:       ## more than one generator at this ref bus
            ## subtract off what is generated by other gens at this bus
            gen[on[refgen[0]], PG] = \
                gen[on[refgen[0]], PG] - sum(gen[on[refgen[1:len(refgen)]], PG])

    ##----- update/compute branch power flows -----
    out = find(branch[:, BR_STATUS] == 0)        ## out-of-service branches
    br =  find(branch[:, BR_STATUS]).astype(int) ## in-service branches

    ## complex power at "from" bus
    Sf = V[ branch[br, F_BUS].astype(int) ] * conj(Yf[br, :] * V) * baseMVA
    ## complex power injected at "to" bus
    St = V[ branch[br, T_BUS].astype(int) ] * conj(Yt[br, :] * V) * baseMVA
    branch[ ix_(br, [PF, QF, PT, QT]) ] = c_[Sf.real, Sf.imag, St.real, St.imag]
    branch[ ix_(out, [PF, QF, PT, QT]) ] = zeros((len(out), 4))

    return bus, gen, branch