numpy.ix_

Here are the examples of the python api numpy.ix_ taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

33 Examples 7

Example 1

Project: attention-lvcsr Source File: test_blocksparse.py
    @staticmethod
    def gemv_numpy2(o, W, h, iIdx, oIdx):
        """
        Other implementation
        """
        from numpy import ix_
        for b in range(o.shape[0]):
            w = W[ix_(iIdx[b], oIdx[b])].swapaxes(1, 2)
            w = w.reshape((w.shape[0] * w.shape[1], w.shape[2] * w.shape[3]))
            o[b] += numpy.dot(h[b].ravel(), w).reshape(o.shape[1:])
        return o

Example 2

Project: attention-lvcsr Source File: test_blocksparse.py
    @staticmethod
    def gemv_numpy3(o, W, h, iIdx, oIdx):
        """
        Other implementation
        """
        from numpy import ix_
        for b in range(o.shape[0]):
            w = W[ix_(iIdx[b], oIdx[b])]
            # The next three lines do the same operation. The last one is the
            # fastest
            # o[b] += (h[b][:, None, :, None] * w).sum(axis=(0, 2))
            # o[b] += numpy.tensordot(h[b], w, [(0,1),(0,2)])
            o[b] += numpy.einsum('ik,ijkl', h[b], w)
        return o

Example 3

Project: pyhawkes Source File: network.py
Function: p
    @property
    def P(self):
        """
        Get the KxK matrix of probabilities
        :return:
        """
        P = self.p[np.ix_(self.c, self.c)]
        if not self.allow_self_connections:
            np.fill_diagonal(P, 0.0)
        return P

Example 4

Project: pyhawkes Source File: network.py
Function: v
    @property
    def V(self):
        """
        Get the KxK matrix of scales
        :return:
        """
        return self.v[np.ix_(self.c, self.c)]

Example 5

Project: pyhawkes Source File: network.py
    def resample_p(self, A):
        """
        Resample p given observations of the weights
        """
        for c1 in xrange(self.C):
            for c2 in xrange(self.C):
                Ac1c2 = A[np.ix_(self.c==c1, self.c==c2)]

                if not self.allow_self_connections:
                    # TODO: Account for self connections
                    pass

                tau1 = self.tau1 + Ac1c2.sum()
                tau0 = self.tau0 + (1-Ac1c2).sum()
                self.p[c1,c2] = np.random.beta(tau1, tau0)

Example 6

Project: pyhawkes Source File: network.py
Function: resample_v
    def resample_v(self, A, W):
        """
        Resample v given observations of the weights
        """
        # import pdb; pdb.set_trace()
        for c1 in xrange(self.C):
            for c2 in xrange(self.C):
                Ac1c2 = A[np.ix_(self.c==c1, self.c==c2)]
                Wc1c2 = W[np.ix_(self.c==c1, self.c==c2)]
                alpha = self.alpha + Ac1c2.sum() * self.kappa
                beta  = self.beta + Wc1c2[Ac1c2 > 0].sum()
                self.v[c1,c2] = np.random.gamma(alpha, 1.0/beta)

Example 7

Project: pyphi Source File: utils.py
def block_reducible(cm, nodes1, nodes2):
    """Return whether connections from ``nodes1`` to ``nodes2`` are reducible.

    Args:
        cm (np.ndarray): The network's connectivity matrix.
        nodes1 (tuple[int]): Source nodes
        nodes2 (tuple[int]): Sink nodes
    """
    if not nodes1 or not nodes2:
        return True  # trivially

    cm = cm[np.ix_(nodes1, nodes2)]

    # Validate the connectivity matrix.
    if not cm.sum(0).all() or not cm.sum(1).all():
        return True
    if len(nodes1) > 1 and len(nodes2) > 1:
        return block_cm(cm)
    return False

Example 8

Project: pyphi Source File: utils.py
def strongly_connected(cm, nodes=None):
    """Return whether the connectivity matrix is strongly connected.

    Args:
        cm (np.ndarray): A square connectivity matrix.

    Keyword Args:
        nodes (tuple[int]): An optional subset of node indices to test strong
            connectivity over.
    """
    if nodes is not None:
        cm = cm[np.ix_(nodes, nodes)]

    num_components, _ = connected_components(cm, connection='strong')
    return num_components < 2

Example 9

Project: holoviews Source File: util.py
Function: cartesian_product
def cartesian_product(arrays):
    """
    Computes the cartesian product of a list of 1D arrays
    returning arrays matching the shape defined by all
    supplied dimensions.
    """
    return np.broadcast_arrays(*np.ix_(*arrays))

Example 10

Project: pgmpy Source File: JointGaussianDistribution.py
Function: marginalize
    def marginalize(self, variables, inplace=True):
        """
        Modifies the distribution with marginalized values.

        Parameters
        ----------

        variables: iterator
                List of variables over which marginalization is to be done.

        inplace: boolean
                If inplace=True it will modify the distribution itself,
                else would return a new distribution.

        Returns
        -------
        JointGaussianDistribution or None :
                if inplace=True (default) returns None
                if inplace=False return a new JointGaussianDistribution instance

        Examples
        --------
        >>> import numpy as np
        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
        >>> dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
        ...             np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
        >>> dis.variables
        ['x1', 'x2', 'x3']
        >>> dis.mean
        array([[ 1],
                [-3],
                [ 4]])
        >>> dis.covariance
        array([[ 4,  2, -2],
               [ 2,  5, -5],
               [-2, -5,  8]])

        >>> dis.marginalize(['x3'])
        dis.variables
        ['x1', 'x2']
        >>> dis.mean
        array([[ 1],
                [-3]]))
        >>> dis.covariance
        narray([[4, 2],
               [2, 5]])
        """

        if not isinstance(variables, list):
            raise TypeError("variables: Expected type list or array-like,\
                             got type {var_type}".format(var_type=type(variables)))

        phi = self if inplace else self.copy()

        var_indexes = [phi.variables.index(var) for var in variables]
        index_to_keep = [self.variables.index(var) for var in self.variables
                         if var not in variables]

        phi.variables = [phi.variables[index] for index in index_to_keep]
        phi.mean = phi.mean[index_to_keep]
        phi.covariance = phi.covariance[np.ix_(index_to_keep, index_to_keep)]
        phi._precision_matrix = None

        if not inplace:
            return phi

Example 11

Project: RLScore Source File: davis_data.py
def setting4_split():
    random.seed(1)
    XD, XT, Y = load_davis()
    drug_ind = range(Y.shape[0])
    target_ind = range(Y.shape[1])
    random.shuffle(drug_ind)
    random.shuffle(target_ind)
    train_drug_ind = drug_ind[:40]
    test_drug_ind = drug_ind[40:]
    train_target_ind = target_ind[:300]
    test_target_ind = target_ind[300:]
    #Setting 4: ensure that d,t pairs do not overlap between
    #training and test set
    Y_train = Y[np.ix_(train_drug_ind, train_target_ind)]
    Y_test = Y[np.ix_(test_drug_ind, test_target_ind)]
    Y_train = Y_train.ravel(order='F')
    Y_test = Y_test.ravel(order='F')
    XD_train = XD[train_drug_ind]
    XT_train = XT[train_target_ind]
    XD_test = XD[test_drug_ind]
    XT_test = XT[test_target_ind]
    return XD_train, XT_train, Y_train, XD_test, XT_test, Y_test    

Example 12

Project: RLScore Source File: global_rankrls.py
    def _reference(self, pairs):
        
        evals, evecs = self.evals, self.svecs
        Y = self.Y
        m = self.size
        
        
        results = []
        
        D = mat(zeros((self.size, 1), dtype=float64))
        C = mat(zeros((self.size, 3), dtype=float64))
        for i in range(self.size):
            D[i, 0] = self.size - 2.
            C[i, 0] = 1.
        lamb = self.regparam
        
        G = evecs * multiply(multiply(evals, 1. / ((m - 2.) * evals + lamb)).T, evecs.T)
        
        
        DY = multiply(D, Y)
        GDY = G * DY
        GC = G * C
        CTG = GC.T
        CT = C.T
        CTGC = CT * GC
        CTY = CT * Y
        CTGDY = CT * GDY
        
        minusI3 = -mat(eye(3))
        for i, j in pairs:
            hoinds = [i, j]
            
            R = mat(zeros((2, 3), dtype=float64))
            R[0, 0] = -1.
            R[1, 0] = -1.
            R[0, 1] = sqrt(self.size - 2.)
            R[1, 2] = sqrt(self.size - 2.)
            
            GBho = GC[hoinds] + G[np.ix_(hoinds, hoinds)] * R
            
            BTGB = CTGC \
                + R.T * GC[hoinds] \
                + CTG[:, hoinds] * R \
                + R.T * G[np.ix_(hoinds, hoinds)] * R
            
            BTY = CTY + R.T * Y[hoinds]
            
            GLYho = GDY[hoinds] - GBho * BTY
            BTGLY = CTGDY + R.T * GDY[hoinds] - BTGB * BTY
            
            F = GLYho - GBho * la.inv(minusI3 + BTGB) * BTGLY
            
            results.append(F)
        return results

Example 13

Project: zipline Source File: frame.py
    def load_adjusted_array(self, columns, dates, assets, mask):
        """
        Load data from our stored baseline.
        """
        column = self.column
        if len(columns) != 1:
            raise ValueError(
                "Can't load multiple columns with DataFrameLoader"
            )
        elif columns[0] != column:
            raise ValueError("Can't load unknown column %s" % columns[0])

        date_indexer = self.dates.get_indexer(dates)
        assets_indexer = self.assets.get_indexer(assets)

        # Boolean arrays with True on matched entries
        good_dates = (date_indexer != -1)
        good_assets = (assets_indexer != -1)

        return {
            column: AdjustedArray(
                # Pull out requested columns/rows from our baseline data.
                data=self.baseline[ix_(date_indexer, assets_indexer)],
                # Mask out requested columns/rows that didnt match.
                mask=(good_assets & as_column(good_dates)) & mask,
                adjustments=self.format_adjustments(dates, assets),
                missing_value=column.missing_value,
            ),
        }

Example 14

Project: qutip Source File: test_sparse.py
def _permutateIndexes(array, row_perm, col_perm):
    return array[np.ix_(row_perm, col_perm)]

Example 15

Project: PYPOWER Source File: makePTDF.py
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

Example 16

Project: PYPOWER Source File: opf.py
def opf(*args):
    """Solves an optimal power flow.

    Returns a C{results} dict.

    The data for the problem can be specified in one of three ways:
      1. a string (ppc) containing the file name of a PYPOWER case
      which defines the data matrices baseMVA, bus, gen, branch, and
      gencost (areas is not used at all, it is only included for
      backward compatibility of the API).
      2. a dict (ppc) containing the data matrices as fields.
      3. the individual data matrices themselves.

    The optional user parameters for user constraints (C{A, l, u}), user costs
    (C{N, fparm, H, Cw}), user variable initializer (C{z0}), and user variable
    limits (C{zl, zu}) can also be specified as fields in a case dict,
    either passed in directly or defined in a case file referenced by name.

    When specified, C{A, l, u} represent additional linear constraints on the
    optimization variables, C{l <= A*[x z] <= u}. If the user specifies an C{A}
    matrix that has more columns than the number of "C{x}" (OPF) variables,
    then there are extra linearly constrained "C{z}" variables. For an
    explanation of the formulation used and instructions for forming the
    C{A} matrix, see the MATPOWER manual.

    A generalized cost on all variables can be applied if input arguments
    C{N}, C{fparm}, C{H} and C{Cw} are specified. First, a linear transformation
    of the optimization variables is defined by means of C{r = N * [x z]}.
    Then, to each element of C{r} a function is applied as encoded in the
    C{fparm} matrix (see MATPOWER manual). If the resulting vector is named
    C{w}, then C{H} and C{Cw} define a quadratic cost on w:
    C{(1/2)*w'*H*w + Cw * w}. C{H} and C{N} should be sparse matrices and C{H}
    should also be symmetric.

    The optional C{ppopt} vector specifies PYPOWER options. If the OPF
    algorithm is not explicitly set in the options PYPOWER will use the default
    solver, based on a primal-dual interior point method. For the AC OPF this
    is C{OPF_ALG = 560}. For the DC OPF, the default is C{OPF_ALG_DC = 200}.
    See L{ppoption} for more details on the available OPF solvers and other OPF
    options and their default values.

    The solved case is returned in a single results dict (described
    below). Also returned are the final objective function value (C{f}) and a
    flag which is C{True} if the algorithm was successful in finding a solution
    (success). Additional optional return values are an algorithm specific
    return status (C{info}), elapsed time in seconds (C{et}), the constraint
    vector (C{g}), the Jacobian matrix (C{jac}), and the vector of variables
    (C{xr}) as well as the constraint multipliers (C{pimul}).

    The single results dict is a PYPOWER case struct (ppc) with the
    usual baseMVA, bus, branch, gen, gencost fields, along with the
    following additional fields:

        - C{order}      see 'help ext2int' for details of this field
        - C{et}         elapsed time in seconds for solving OPF
        - C{success}    1 if solver converged successfully, 0 otherwise
        - C{om}         OPF model object, see 'help opf_model'
        - C{x}          final value of optimization variables (internal order)
        - C{f}          final objective function value
        - C{mu}         shadow prices on ...
            - C{var}
                - C{l}  lower bounds on variables
                - C{u}  upper bounds on variables
            - C{nln}
                - C{l}  lower bounds on nonlinear constraints
                - C{u}  upper bounds on nonlinear constraints
            - C{lin}
                - C{l}  lower bounds on linear constraints
                - C{u}  upper bounds on linear constraints
        - C{g}          (optional) constraint values
        - C{dg}         (optional) constraint 1st derivatives
        - C{df}         (optional) obj fun 1st derivatives (not yet implemented)
        - C{d2f}        (optional) obj fun 2nd derivatives (not yet implemented)
        - C{raw}        raw solver output in form returned by MINOS, and more
            - C{xr}     final value of optimization variables
            - C{pimul}  constraint multipliers
            - C{info}   solver specific termination code
            - C{output} solver specific output information
               - C{alg} algorithm code of solver used
        - C{var}
            - C{val}    optimization variable values, by named block
                - C{Va}     voltage angles
                - C{Vm}     voltage magnitudes (AC only)
                - C{Pg}     real power injections
                - C{Qg}     reactive power injections (AC only)
                - C{y}      constrained cost variable (only if have pwl costs)
                - (other) any user defined variable blocks
            - C{mu}     variable bound shadow prices, by named block
                - C{l}  lower bound shadow prices
                    - C{Va}, C{Vm}, C{Pg}, C{Qg}, C{y}, (other)
                - C{u}  upper bound shadow prices
                    - C{Va}, C{Vm}, C{Pg}, C{Qg}, C{y}, (other)
        - C{nln}    (AC only)
            - C{mu}     shadow prices on nonlinear constraints, by named block
                - C{l}  lower bounds
                    - C{Pmis}   real power mismatch equations
                    - C{Qmis}   reactive power mismatch equations
                    - C{Sf}     flow limits at "from" end of branches
                    - C{St}     flow limits at "to" end of branches
                - C{u}  upper bounds
                    - C{Pmis}, C{Qmis}, C{Sf}, C{St}
        - C{lin}
            - C{mu}     shadow prices on linear constraints, by named block
                - C{l}  lower bounds
                    - C{Pmis}   real power mistmatch equations (DC only)
                    - C{Pf}     flow limits at "from" end of branches (DC only)
                    - C{Pt}     flow limits at "to" end of branches (DC only)
                    - C{PQh}    upper portion of gen PQ-capability curve(AC only)
                    - C{PQl}    lower portion of gen PQ-capability curve(AC only)
                    - C{vl}     constant power factor constraint for loads
                    - C{ycon}   basin constraints for CCV for pwl costs
                    - (other) any user defined constraint blocks
                - C{u}  upper bounds
                    - C{Pmis}, C{Pf}, C{Pf}, C{PQh}, C{PQl}, C{vl}, C{ycon},
                    - (other)
        - C{cost}       user defined cost values, by named block

    @see: L{runopf}, L{dcopf}, L{uopf}, L{caseformat}

    @author: Ray Zimmerman (PSERC Cornell)
    @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad
    Autonoma de Manizales)
    """
    ##----- initialization -----
    t0 = time()         ## start timer

    ## process input arguments
    ppc, ppopt = opf_args2(*args)

    ## add zero columns to bus, gen, branch for multipliers, etc if needed
    nb   = shape(ppc['bus'])[0]    ## number of buses
    nl   = shape(ppc['branch'])[0] ## number of branches
    ng   = shape(ppc['gen'])[0]    ## number of dispatchable injections
    if shape(ppc['bus'])[1] < MU_VMIN + 1:
        ppc['bus'] = c_[ppc['bus'], zeros((nb, MU_VMIN + 1 - shape(ppc['bus'])[1]))]

    if shape(ppc['gen'])[1] < MU_QMIN + 1:
        ppc['gen'] = c_[ppc['gen'], zeros((ng, MU_QMIN + 1 - shape(ppc['gen'])[1]))]

    if shape(ppc['branch'])[1] < MU_ANGMAX + 1:
        ppc['branch'] = c_[ppc['branch'], zeros((nl, MU_ANGMAX + 1 - shape(ppc['branch'])[1]))]

    ##-----  convert to internal numbering, remove out-of-service stuff  -----
    ppc = ext2int(ppc)

    ##-----  construct OPF model object  -----
    om = opf_setup(ppc, ppopt)

    ##-----  execute the OPF  -----
    results, success, raw = opf_execute(om, ppopt)

    ##-----  revert to original ordering, including out-of-service stuff  -----
    results = int2ext(results)

    ## zero out result fields of out-of-service gens & branches
    if len(results['order']['gen']['status']['off']) > 0:
        results['gen'][ ix_(results['order']['gen']['status']['off'], [PG, QG, MU_PMAX, MU_PMIN]) ] = 0

    if len(results['order']['branch']['status']['off']) > 0:
        results['branch'][ ix_(results['order']['branch']['status']['off'], [PF, QF, PT, QT, MU_SF, MU_ST, MU_ANGMIN, MU_ANGMAX]) ] = 0

    ##-----  finish preparing output  -----
    et = time() - t0      ## compute elapsed time

    results['et'] = et
    results['success'] = success
    results['raw'] = raw

    return results

Example 17

Project: PYPOWER Source File: pfsoln.py
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)
    """
    ## 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

Example 18

Project: PYPOWER Source File: runpf.py
def runpf(casedata=None, ppopt=None, fname='', solvedcase=''):
    """Runs a power flow.

    Runs a power flow [full AC Newton's method by default] and optionally
    returns the solved values in the data matrices, a flag which is C{True} if
    the algorithm was successful in finding a solution, and the elapsed
    time in seconds. All input arguments are optional. If C{casename} is
    provided it specifies the name of the input data file or dict
    containing the power flow data. The default value is 'case9'.

    If the ppopt is provided it overrides the default PYPOWER options
    vector and can be used to specify the solution algorithm and output
    options among other things. If the 3rd argument is given the pretty
    printed output will be appended to the file whose name is given in
    C{fname}. If C{solvedcase} is specified the solved case will be written
    to a case file in PYPOWER format with the specified name. If C{solvedcase}
    ends with '.mat' it saves the case as a MAT-file otherwise it saves it
    as a Python-file.

    If the C{ENFORCE_Q_LIMS} options is set to C{True} [default is false] then
    if any generator reactive power limit is violated after running the AC
    power flow, the corresponding bus is converted to a PQ bus, with Qg at
    the limit, and the case is re-run. The voltage magnitude at the bus
    will deviate from the specified value in order to satisfy the reactive
    power limit. If the reference bus is converted to PQ, the first
    remaining PV bus will be used as the slack bus for the next iteration.
    This may result in the real power output at this generator being
    slightly off from the specified values.

    Enforcing of generator Q limits inspired by contributions from Mu Lin,
    Lincoln University, New Zealand (1/14/05).

    @author: Ray Zimmerman (PSERC Cornell)
    """
    ## default arguments
    if casedata is None:
        casedata = join(dirname(__file__), 'case9')
    ppopt = ppoption(ppopt)

    ## options
    verbose = ppopt["VERBOSE"]
    qlim = ppopt["ENFORCE_Q_LIMS"]  ## enforce Q limits on gens?
    dc = ppopt["PF_DC"]             ## use DC formulation?

    ## read data
    ppc = loadcase(casedata)

    ## add zero columns to branch for flows if needed
    if ppc["branch"].shape[1] < QT:
        ppc["branch"] = c_[ppc["branch"],
                           zeros((ppc["branch"].shape[0],
                                  QT - ppc["branch"].shape[1] + 1))]

    ## convert to internal indexing
    ppc = ext2int(ppc)
    baseMVA, bus, gen, branch = \
        ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"]

    ## get bus index lists of each type of bus
    ref, pv, pq = bustypes(bus, gen)

    ## generator info
    on = find(gen[:, GEN_STATUS] > 0)      ## which generators are on?
    gbus = gen[on, GEN_BUS].astype(int)    ## what buses are they at?

    ##-----  run the power flow  -----
    t0 = time()
    if verbose > 0:
        v = ppver('all')
        stdout.write('PYPOWER Version %s, %s' % (v["Version"], v["Date"]))

    if dc:                               # DC formulation
        if verbose:
            stdout.write(' -- DC Power Flow\n')

        ## initial state
        Va0 = bus[:, VA] * (pi / 180)

        ## build B matrices and phase shift injections
        B, Bf, Pbusinj, Pfinj = makeBdc(baseMVA, bus, branch)

        ## compute complex bus power injections [generation - load]
        ## adjusted for phase shifters and real shunts
        Pbus = makeSbus(baseMVA, bus, gen).real - Pbusinj - bus[:, GS] / baseMVA

        ## "run" the power flow
        Va = dcpf(B, Pbus, Va0, ref, pv, pq)

        ## update data matrices with solution
        branch[:, [QF, QT]] = zeros((branch.shape[0], 2))
        branch[:, PF] = (Bf * Va + Pfinj) * baseMVA
        branch[:, PT] = -branch[:, PF]
        bus[:, VM] = ones(bus.shape[0])
        bus[:, VA] = Va * (180 / pi)
        ## update Pg for slack generator (1st gen at ref bus)
        ## (note: other gens at ref bus are accounted for in Pbus)
        ##      Pg = Pinj + Pload + Gs
        ##      newPg = oldPg + newPinj - oldPinj
        refgen = zeros(len(ref), dtype=int)
        for k in range(len(ref)):
            temp = find(gbus == ref[k])
            refgen[k] = on[temp[0]]
        gen[refgen, PG] = gen[refgen, PG] + (B[ref, :] * Va - Pbus[ref]) * baseMVA

        success = 1
    else:                                ## AC formulation
        alg = ppopt['PF_ALG']
        if verbose > 0:
            if alg == 1:
                solver = 'Newton'
            elif alg == 2:
                solver = 'fast-decoupled, XB'
            elif alg == 3:
                solver = 'fast-decoupled, BX'
            elif alg == 4:
                solver = 'Gauss-Seidel'
            else:
                solver = 'unknown'
            print(' -- AC Power Flow (%s)\n' % solver)

        ## initial state
        # V0    = ones(bus.shape[0])            ## flat start
        V0  = bus[:, VM] * exp(1j * pi/180 * bus[:, VA])
        V0[gbus] = gen[on, VG] / abs(V0[gbus]) * V0[gbus]

        if qlim:
            ref0 = ref                         ## save index and angle of
            Varef0 = bus[ref0, VA]             ##   original reference bus(es)
            limited = []                       ## list of indices of gens @ Q lims
            fixedQg = zeros(gen.shape[0])      ## Qg of gens at Q limits

        repeat = True
        while repeat:
            ## build admittance matrices
            Ybus, Yf, Yt = makeYbus(baseMVA, bus, branch)

            ## compute complex bus power injections [generation - load]
            Sbus = makeSbus(baseMVA, bus, gen)

            ## run the power flow
            alg = ppopt["PF_ALG"]
            if alg == 1:
                V, success, _ = newtonpf(Ybus, Sbus, V0, ref, pv, pq, ppopt)
            elif alg == 2 or alg == 3:
                Bp, Bpp = makeB(baseMVA, bus, branch, alg)
                V, success, _ = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, ppopt)
            elif alg == 4:
                V, success, _ = gausspf(Ybus, Sbus, V0, ref, pv, pq, ppopt)
            else:
                stderr.write('Only Newton''s method, fast-decoupled, and '
                             'Gauss-Seidel power flow algorithms currently '
                             'implemented.\n')

            ## update data matrices with solution
            bus, gen, branch = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq)

            if qlim:             ## enforce generator Q limits
                ## find gens with violated Q constraints
                gen_status = gen[:, GEN_STATUS] > 0
                qg_max_lim = gen[:, QG] > gen[:, QMAX]
                qg_min_lim = gen[:, QG] < gen[:, QMIN]
                
                mx = find( gen_status & qg_max_lim )
                mn = find( gen_status & qg_min_lim )
                
                if len(mx) > 0 or len(mn) > 0:  ## we have some Q limit violations
                    # No PV generators
                    if len(pv) == 0:
                        if verbose:
                            if len(mx) > 0:
                                print('Gen %d [only one left] exceeds upper Q limit : INFEASIBLE PROBLEM\n' % mx + 1)
                            else:
                                print('Gen %d [only one left] exceeds lower Q limit : INFEASIBLE PROBLEM\n' % mn + 1)

                        success = 0
                        break

                    ## one at a time?
                    if qlim == 2:    ## fix largest violation, ignore the rest
                        k = argmax(r_[gen[mx, QG] - gen[mx, QMAX],
                                      gen[mn, QMIN] - gen[mn, QG]])
                        if k > len(mx):
                            mn = mn[k - len(mx)]
                            mx = []
                        else:
                            mx = mx[k]
                            mn = []

                    if verbose and len(mx) > 0:
                        for i in range(len(mx)):
                            print('Gen ' + str(mx[i] + 1) + ' at upper Q limit, converting to PQ bus\n')

                    if verbose and len(mn) > 0:
                        for i in range(len(mn)):
                            print('Gen ' + str(mn[i] + 1) + ' at lower Q limit, converting to PQ bus\n')

                    ## save corresponding limit values
                    fixedQg[mx] = gen[mx, QMAX]
                    fixedQg[mn] = gen[mn, QMIN]
                    mx = r_[mx, mn].astype(int)

                    ## convert to PQ bus
                    gen[mx, QG] = fixedQg[mx]      ## set Qg to binding 
                    for i in range(len(mx)):            ## [one at a time, since they may be at same bus]
                        gen[mx[i], GEN_STATUS] = 0        ## temporarily turn off gen,
                        bi = gen[mx[i], GEN_BUS]   ## adjust load accordingly,
                        bus[bi, [PD, QD]] = (bus[bi, [PD, QD]] - gen[mx[i], [PG, QG]])
                    
                    if len(ref) > 1 and any(bus[gen[mx, GEN_BUS], BUS_TYPE] == REF):
                        raise ValueError('Sorry, PYPOWER cannot enforce Q '
                                         'limits for slack buses in systems '
                                         'with multiple slacks.')
                    
                    bus[gen[mx, GEN_BUS].astype(int), BUS_TYPE] = PQ   ## & set bus type to PQ

                    ## update bus index lists of each type of bus
                    ref_temp = ref
                    ref, pv, pq = bustypes(bus, gen)
                    if verbose and ref != ref_temp:
                        print('Bus %d is new slack bus\n' % ref)

                    limited = r_[limited, mx].astype(int)
                else:
                    repeat = 0 ## no more generator Q limits violated
            else:
                repeat = 0     ## don't enforce generator Q limits, once is enough

        if qlim and len(limited) > 0:
            ## restore injections from limited gens [those at Q limits]
            gen[limited, QG] = fixedQg[limited]    ## restore Qg value,
            for i in range(len(limited)):               ## [one at a time, since they may be at same bus]
                bi = gen[limited[i], GEN_BUS]           ## re-adjust load,
                bus[bi, [PD, QD]] = bus[bi, [PD, QD]] + gen[limited[i], [PG, QG]]
                gen[limited[i], GEN_STATUS] = 1           ## and turn gen back on
            
            if ref != ref0:
                ## adjust voltage angles to make original ref bus correct
                bus[:, VA] = bus[:, VA] - bus[ref0, VA] + Varef0

    ppc["et"] = time() - t0
    ppc["success"] = success

    ##-----  output results  -----
    ## convert back to original bus numbering & print results
    ppc["bus"], ppc["gen"], ppc["branch"] = bus, gen, branch
    results = int2ext(ppc)

    ## zero out result fields of out-of-service gens & branches
    if len(results["order"]["gen"]["status"]["off"]) > 0:
        results["gen"][ix_(results["order"]["gen"]["status"]["off"], [PG, QG])] = 0

    if len(results["order"]["branch"]["status"]["off"]) > 0:
        results["branch"][ix_(results["order"]["branch"]["status"]["off"], [PF, QF, PT, QT])] = 0

    if fname:
        fd = None
        try:
            fd = open(fname, "a")
        except Exception as detail:
            stderr.write("Error opening %s: %s.\n" % (fname, detail))
        finally:
            if fd is not None:
                printpf(results, fd, ppopt)
                fd.close()
    else:
        printpf(results, stdout, ppopt)

    ## save solved case
    if solvedcase:
        savecase(solvedcase, results)

    return results, success

Example 19

Project: PYPOWER Source File: scale_load.py
def scale_load(load, bus, gen=None, load_zone=None, opt=None):
    """Scales fixed and/or dispatchable loads.

    Assumes consecutive bus numbering when dealing with dispatchable loads.

    @param load: Each element specifies the amount of scaling for the
        corresponding load zone, either as a direct scale factor
        or as a target quantity. If there are C{nz} load zones this
        vector has C{nz} elements.
    @param bus: Standard C{bus} matrix with C{nb} rows, where the fixed active
        and reactive loads available for scaling are specified in
        columns C{PD} and C{QD}
    @param gen: (optional) standard C{gen} matrix with C{ng} rows, where the
        dispatchable loads available for scaling are specified by
        columns C{PG}, C{QG}, C{PMIN}, C{QMIN} and C{QMAX} (in rows for which
        C{isload(gen)} returns C{true}). If C{gen} is empty, it assumes
        there are no dispatchable loads.
    @param load_zone: (optional) C{nb} element vector where the value of
        each element is either zero or the index of the load zone
        to which the corresponding bus belongs. If C{load_zone[b] = k}
        then the loads at bus C{b} will be scaled according to the
        value of C{load[k]}. If C{load_zone[b] = 0}, the loads at bus C{b}
        will not be modified. If C{load_zone} is empty, the default is
        determined by the dimensions of the C{load} vector. If C{load} is
        a scalar, a single system-wide zone including all buses is
        used, i.e. C{load_zone = ones(nb)}. If C{load} is a vector, the
        default C{load_zone} is defined as the areas specified in the
        C{bus} matrix, i.e. C{load_zone = bus[:, BUS_AREA]}, and C{load}
        should have dimension C{= max(bus[:, BUS_AREA])}.
    @param opt: (optional) dict with three possible fields, 'scale',
        'pq' and 'which' that determine the behavior as follows:
            - C{scale} (default is 'FACTOR')
                - 'FACTOR'   : C{load} consists of direct scale factors, where
                C{load[k] =} scale factor C{R[k]} for zone C{k}
                - 'QUANTITY' : C{load} consists of target quantities, where
                C{load[k] =} desired total active load in MW for
                zone C{k} after scaling by an appropriate C{R(k)}
            - C{pq}    (default is 'PQ')
                - 'PQ' : scale both active and reactive loads
                - 'P'  : scale only active loads
            - C{which} (default is 'BOTH' if GEN is provided, else 'FIXED')
                - 'FIXED'        : scale only fixed loads
                - 'DISPATCHABLE' : scale only dispatchable loads
                - 'BOTH'         : scale both fixed and dispatchable loads

    @see: L{total_load}

    @author: Ray Zimmerman (PSERC Cornell)
    """
    nb = bus.shape[0]   ## number of buses

    ##-----  process inputs  -----
    bus = bus.copy()
    if gen is None:
        gen = array([])
    else:
        gen = gen.copy()
    if load_zone is None:
        load_zone = array([], int)
    if opt is None:
        opt = {}

    ## fill out and check opt
    if len(gen) == 0:
        opt["which"] = 'FIXED'
    if 'pq' not in opt:
        opt["pq"] = 'PQ'          ## 'PQ' or 'P'
    if 'which' not in opt:
        opt["which"] = 'BOTH'     ## 'FIXED', 'DISPATCHABLE' or 'BOTH'
    if 'scale' not in opt:
        opt["scale"] = 'FACTOR'   ## 'FACTOR' or 'QUANTITY'
    if (opt["pq"] != 'P') and (opt["pq"] != 'PQ'):
        stderr.write("scale_load: opt['pq'] must equal 'PQ' or 'P'\n")
    if (opt["which"][0] != 'F') and (opt["which"][0] != 'D') and (opt["which"][0] != 'B'):
        stderr.write("scale_load: opt.which should be 'FIXED, 'DISPATCHABLE or 'BOTH'\n")
    if (opt["scale"][0] != 'F') and (opt["scale"][0] != 'Q'):
        stderr.write("scale_load: opt.scale should be 'FACTOR or 'QUANTITY'\n")
    if (len(gen) == 0) and (opt["which"][0] != 'F'):
        stderr.write('scale_load: need gen matrix to scale dispatchable loads\n')

    ## create dispatchable load connection matrix
    if len(gen) > 0:
        ng = gen.shape[0]
        is_ld = isload(gen) & (gen[:, GEN_STATUS] > 0)
        ld = find(is_ld)

        ## create map of external bus numbers to bus indices
        i2e = bus[:, BUS_I].astype(int)
        e2i = zeros(max(i2e) + 1, int)
        e2i[i2e] = arange(nb)

        gbus = gen[:, GEN_BUS].astype(int)
        Cld = sparse((is_ld, (e2i[gbus], arange(ng))), (nb, ng))
    else:
        ng = 0
        ld = array([], int)

    if len(load_zone) == 0:
        if len(load) == 1:        ## make a single zone of all load buses
            load_zone = zeros(nb, int)             ## initialize
            load_zone[bus[:, PD] != 0 or bus[:, QD] != 0] = 1  ## FIXED loads
            if len(gen) > 0:
                gbus = gen[ld, GEN_BUS].astype(int)
                load_zone[e2i[gbus]] = 1    ## DISPATCHABLE loads
        else:                        ## use areas defined in bus data as zones
            load_zone = bus[:, BUS_AREA]

    ## check load_zone to make sure it's consistent with size of load vector
    if max(load_zone) > len(load):
        stderr.write('scale_load: load vector must have a value for each load zone specified\n')

    ##-----  compute scale factors for each zone  -----
    scale = load.copy()
    Pdd = zeros(nb)     ## dispatchable P at each bus
    if opt["scale"][0] == 'Q':  ## 'QUANTITY'
        ## find load capacity from dispatchable loads
        if len(gen) > 0:
            Pdd = -Cld * gen[:, PMIN]

        ## compute scale factors
        for k in range(len(load)):
            idx = find(load_zone == k + 1)
            fixed = sum(bus[idx, PD])
            dispatchable = sum(Pdd[idx])
            total = fixed + dispatchable
            if opt["which"][0] == 'B':      ## 'BOTH'
                if total != 0:
                    scale[k] = load[k] / total
                elif load[k] == total:
                    scale[k] = 1
                else:
                    raise ScalingError('scale_load: impossible to make zone %d load equal %g by scaling non-existent loads\n' % (k, load[k]))
            elif opt["which"][0] == 'F':    ## 'FIXED'
                if fixed != 0:
                    scale[k] = (load[k] - dispatchable) / fixed
                elif load[k] == dispatchable:
                    scale[k] = 1
                else:
                    raise ScalingError('scale_load: impossible to make zone %d load equal %g by scaling non-existent fixed load\n' % (k, load[k]))
            elif opt["which"][0] == 'D':    ## 'DISPATCHABLE'
                if dispatchable != 0:
                    scale[k] = (load[k] - fixed) / dispatchable
                elif load[k] == fixed:
                    scale[k] = 1
                else:
                    raise ScalingError('scale_load: impossible to make zone %d load equal %g by scaling non-existent dispatchable load\n' % (k, load[k]))

    ##-----  do the scaling  -----
    ## fixed loads
    if opt["which"][0] != 'D':      ## includes 'FIXED', not 'DISPATCHABLE' only
        for k in range(len(scale)):
            idx = find(load_zone == k + 1)
            bus[idx, PD] = bus[idx, PD] * scale[k]
            if opt["pq"] == 'PQ':
                bus[idx, QD] = bus[idx, QD] * scale[k]

    ## dispatchable loads
    if opt["which"][0] != 'F':      ## includes 'DISPATCHABLE', not 'FIXED' only
        for k in range(len(scale)):
            idx = find(load_zone == k + 1)
            gbus = gen[ld, GEN_BUS].astype(int)
            i = find( in1d(e2i[gbus], idx) )
            ig = ld[i]

            gen[ix_(ig, [PG, PMIN])] = gen[ix_(ig, [PG, PMIN])] * scale[k]
            if opt["pq"] == 'PQ':
                gen[ix_(ig, [QG, QMIN, QMAX])] = gen[ix_(ig, [QG, QMIN, QMAX])] * scale[k]

    return bus, gen

Example 20

Project: PYPOWER Source File: t_off2case.py
def t_off2case(quiet=False):
    """Tests for code in C{off2case}.

    @author: Ray Zimmerman (PSERC Cornell)
    """
    n_tests = 35

    t_begin(n_tests, quiet)

    ## generator data
    #    bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
    gen0 = array([
        [1,   10,   0,  60, -15, 1, 100, 1, 60, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [2,   10,   0,  60, -15, 1, 100, 1, 60, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [7,  -30, -15,   0, -15, 1, 100, 1, 0, -30, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [13,  10,   0,  60, -15, 1, 100, 1, 60, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [30, -30, 7.5, 7.5,   0, 1, 100, 1, 0, -30, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    ], float)
    ## generator cost data
    #    1    startup    shutdown    n    x1    y1    ...    xn    yn
    #    2    startup    shutdown    n    c(n-1)    ...    c0
    gencost0 = array([
        [1, 0,   0, 4,   0, 0,  12, 240,   36, 1200, 60, 2400],
        [1, 100, 0, 4,   0, 0,  12, 240,   36, 1200, 60, 2400],
        [1, 0,   0, 4, -30, 0, -20, 1000, -10, 2000,  0, 3000],
        [1, 0,   0, 4,   0, 0,  12, 240,   36, 1200, 60, 2400],
        [1, 0,  50, 4, -30, 0, -20, 1000, -10, 2000,  0, 3000]
    ], float)

    try:
        from pypower.extras.smartmarket import off2case
    except ImportError:
        t_skip(n_tests, 'smartmarket code not available')
        return

    t = 'isload()'
    t_is(isload(gen0), array([0, 0, 1, 0, 1], bool), 8, t)

    G = find( ~isload(gen0) )
    L = find(  isload(gen0) )
    nGL = len(G) + len(L)

    t = 'P offers only';
    offers = {'P': {}}
    offers['P']['qty'] = array([[25], [26], [27]], float)
    offers['P']['prc'] = array([[10], [50], [100]], float)
    gen, gencost = off2case(gen0, gencost0, offers)

    gen1 = gen0.copy()
    gen1[G, PMAX] = offers['P']['qty'].flatten()
    gen1[L, GEN_STATUS] = 0
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0.copy()
    gencost1[ix_(G, range(NCOST, NCOST + 9))] = c_[array([
        [2, 0, 0, 25,  250],
        [2, 0, 0, 26, 1300],
        [2, 0, 0, 27, 2700],
    ]), zeros((3, 4))]

    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    offers['P']['qty'] = array([[25], [26], [0], [27],  [0]], float)
    offers['P']['prc'] = array([[10], [50], [0], [100], [0]], float)
    gen, gencost = off2case(gen0, gencost0, offers)
    t_is( gen, gen1, 8, [t, ' (all rows in offer) - gen'] )
    t_is( gencost, gencost1, 8, [t, ' (all rows in offer) - gencost'] )

    t = 'P offers only (GEN_STATUS=0 for 0 qty offer)';
    offers['P']['qty'] = array([ [0], [26],  [27]], float)
    offers['P']['prc'] = array([[10], [50], [100]], float)
    gen, gencost = off2case(gen0, gencost0, offers)

    gen1 = gen0.copy()
    gen1[G[1:3], PMAX] = offers['P']['qty'].flatten()[1:3]
    gen1[G[0], GEN_STATUS] = 0
    gen1[L, GEN_STATUS] = 0
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0.copy()
    gencost1[ix_(G[1:3], range(NCOST, NCOST + 9))] = c_[array([
        [2, 0, 0, 26, 1300],
        [2, 0, 0, 27, 2700]
    ]), zeros((2, 4))]

    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'P offers, lim[\'P\'][\'max_offer\']';
    offers['P']['qty'] = array([[25], [26], [27]], float)
    offers['P']['prc'] = array([[10], [50], [100]], float)
    lim = {'P': {'max_offer': 75}}
    gen, gencost = off2case(gen0, gencost0, offers, lim=lim)

    gen1 = gen0.copy()
    gen1[G[:2], PMAX] = offers['P']['qty'].flatten()[:2, :]
    gen1[r_[G[2], L], GEN_STATUS] = 0
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0.copy()
    gencost1[ix_(G[:2], range(NCOST, NCOST + 9))] = c_[array([
        [2, 0, 0, 25,  250],
        [2, 0, 0, 26, 1300]
    ]), zeros((2, 4))]
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'P offers & P bids';
    bids = {'P': {'qty': array([ [20], [28]], float),
                  'prc': array([[100], [10]], float)}}
    gen, gencost = off2case(gen0, gencost0, offers, bids)

    gen1 = gen0.copy()
    gen1[G, PMAX] = offers['P']['qty']
    gen1[ix_(L, [PMIN, QMIN, QMAX])] = array([
        [-20, -10, 0],
        [-28,   0, 7]
    ])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0[:, :8].copy()
    gencost1[ix_(G, range(NCOST, NCOST + 4))] = array([
        [2, 0, 0, 25,  250],
        [2, 0, 0, 26, 1300],
        [2, 0, 0, 27, 2700]
    ])
    gencost1[ix_(L, range(NCOST, NCOST + 4))] = array([
        [2, -20, -2000, 0, 0],
        [2, -28,  -280, 0, 0]
    ])
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'P offers & P bids (all rows in bid)';
    bids['P']['qty'] = array([[0], [0],  [20], [0], [28]], float)
    bids['P']['prc'] = array([[0], [0], [100], [0], [10]], float)
    gen, gencost = off2case(gen0, gencost0, offers, bids)

    t_is( gen, gen1, 8, [t, ' - gen'] )
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'P offers & P bids (GEN_STATUS=0 for 0 qty bid)';
    bids['P']['qty'] = array([  [0], [28]], float)
    bids['P']['prc'] = array([[100], [10]], float)
    gen, gencost = off2case(gen0, gencost0, offers, bids)

    gen1 = gen0.copy()
    gen1[G, PMAX] = offers['P']['qty']
    gen1[L[0], GEN_STATUS] = 0
    gen1[L[1], [PMIN, QMIN, QMAX]] = array([-28, 0, 7])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0.copy()
    gencost1[ix_(G, range(NCOST, NCOST + 9))] = c_[array([
        [2, 0, 0, 25, 250],
        [2, 0, 0, 26, 1300],
        [2, 0, 0, 27, 2700]
    ]), zeros((3, 4))]
    gencost1[L[1], NCOST:NCOST + 8] = c_[array([
        [2, -28, -280, 0, 0]
    ]), zeros((1, 4))]
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'P offers & P bids (1 gen with both)';
    gen2 = gen0.copy()
    gen2[1, PMIN] = -5
    bids['P']['qty'] = array([[0],  [3],  [20], [0], [28]], float)
    bids['P']['prc'] = array([[0], [50], [100], [0], [10]], float)
    gen, gencost = off2case(gen2, gencost0, offers, bids)

    gen1 = gen2.copy()
    gen1[G, PMAX] = offers['P']['qty']
    gen1[1, PMIN] = -sum( bids['P']['qty'][1, :] )
    gen1[ix_(L, [PMIN, QMIN, QMAX])] = array([
        [-20, -10, 0],
        [-28,   0, 7]
    ])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0[:, :10].copy()
    gencost1[ix_(G, range(NCOST, NCOST + 7))] = array([
        [2,  0,    0, 25,  250,  0,    0],
        [3, -3, -150,  0,    0, 26, 1300],
        [2,  0,    0, 27, 2700,  0,    0]
    ])
    gencost1[ix_(L, range(NCOST, NCOST + 7))] = c_[array([
        [2, -20, -2000, 0, 0],
        [2, -28,  -280, 0, 0]
    ]), zeros((2, 2))]
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'P offers & P bids, lim[\'P\'][\'max_offer\']/[\'min_bid\']'
    bids['P']['qty'] = array([[20],  [28]], float)
    bids['P']['prc'] = array([[100], [10]], float)
    lim['P']['min_bid'] = 50.0
    gen, gencost = off2case(gen0, gencost0, offers, bids, lim)

    gen1 = gen0.copy()
    gen1[G[:2], PMAX] = offers['P']['qty'][:2, :]
    gen1[r_[G[2], L[1]], GEN_STATUS] = 0
    gen1[L[0], [PMIN, QMIN, QMAX]] = array([-20, -10, 0])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0.copy()
    gencost1[ix_(G[:2], range(NCOST, NCOST + 9))] = c_[array([
        [2, 0, 0, 25,  250],
        [2, 0, 0, 26, 1300]
    ]), zeros((2, 4))]
    gencost1[L[0], NCOST:NCOST + 9] = array([2, -20, -2000, 0, 0, 0, 0, 0, 0])
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'P offers & P bids, lim[\'P\'][\'max_offer\']/[\'min_bid\'], multi-block'
    offers['P']['qty'] = array([[10,  40], [20, 30], [25, 25]], float)
    offers['P']['prc'] = array([[10, 100], [25, 65], [50, 90]], float)
    bids['P']['qty'] = array([[ 20, 10], [12, 18]], float)
    bids['P']['prc'] = array([[100, 60], [70, 10]], float)
    gen, gencost = off2case(gen0, gencost0, offers, bids, lim)

    gen1 = gen0.copy()
    gen1[G, PMAX] = array([10, 50, 25])
    gen1[ix_(L, [PMIN, QMIN, QMAX])] = array([
        [-30, -15, 0],
        [-12,   0, 3]
    ])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0[:, :10].copy()
    gencost1[ix_(G, range(NCOST, NCOST + 7))] = array([
        [2, 0, 0, 10,  100, 0,     0],
        [3, 0, 0, 20,  500, 50, 2450],
        [2, 0, 0, 25, 1250, 0,     0]
    ])
    gencost1[ix_(L, range(NCOST, NCOST + 7))] = array([
        [3, -30, -2600, -20, -2000, 0, 0],
        [2, -12,  -840,   0,     0, 0, 0]
    ])
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    ##-----  reactive  -----
    ## generator cost data
    #    1    startup    shutdown    n    x1    y1    ...    xn    yn
    #    2    startup    shutdown    n    c(n-1)    ...    c0
    gencost0 = array([
        [1,   0,  0, 4,   0,    0,  12,  240,  36, 1200, 60, 2400],
        [1, 100,  0, 4,   0,    0,  12,  240,  36, 1200, 60, 2400],
        [1,   0,  0, 4, -30,    0, -20, 1000, -10, 2000,  0, 3000],
        [1,   0,  0, 4,   0,    0,  12,  240,  36, 1200, 60, 2400],
        [1,   0, 50, 4, -30,    0, -20, 1000, -10, 2000,  0, 3000],
        [1,   0,  0, 4, -15, -150,   0,    0,  30,  150, 60,  450],
        [1, 100,  0, 2,   0,    0,   0,    0,   0,    0,  0,    0],
        [1,   0,  0, 3, -20,  -15, -10,  -10,   0,    0,  0,    0],
        [1,   0,  0, 3,   0,    0,  40,   80,  60,  180,  0,    0],
        [1,   0, 50, 2,   0,    0,   0,    0,   0,    0,  0,    0]
    ], float)

    t = 'PQ offers only';
    offers['P']['qty'] = array([[25], [26],  [27]], float)
    offers['P']['prc'] = array([[10], [50], [100]], float)
    offers['Q']['qty'] = array([[10], [20],  [30]], float)
    offers['Q']['prc'] = array([[10],  [5],   [1]], float)
    gen, gencost = off2case(gen0, gencost0, offers)

    gen1 = gen0.copy()
    gen1[G, PMAX] = offers['P']['qty']
    gen1[G, QMAX] = offers['Q']['qty']
    gen1[G, QMIN] = 0
    gen1[L, GEN_STATUS] = 0
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0.copy()
    gencost1[ix_(G, range(NCOST, NCOST + 9))] = c_[array([
        [2, 0, 0, 25,  250],
        [2, 0, 0, 26, 1300],
        [2, 0, 0, 27, 2700]
    ]), zeros((3, 4))]
    gencost1[ix_(G + nGL - 1, range(NCOST, NCOST + 9))] = c_[array([
        [2, 0, 0, 10, 100],
        [2, 0, 0, 20, 100],
        [2, 0, 0, 30,  30]
    ]), zeros((3, 4))]

    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'PQ offers & PQ bids, lim.P/Q.max_offer/min_bid, multi-block';
    offers['P']['qty'] = array([[10,  40], [20, 30], [25, 25]], float)
    offers['P']['prc'] = array([[10, 100], [25, 65], [50, 90]], float)
    bids['P']['qty'] = array([[ 20, 10], [12, 18]], float)
    bids['P']['prc'] = array([[100, 60], [70, 10]], float)
    offers['Q']['qty'] = array([[ 5,  5], [10, 10], [15, 15]], float)
    offers['Q']['prc'] = array([[10, 20], [ 5, 60], [ 1, 10]], float)
    bids['Q']['qty'] = array([ 15, 10, 15,  15,  0], float)
    bids['Q']['prc'] = array([-10,  0,  5, -20, 10], float)
    lim['Q']['max_offer'] = 50.0
    lim['Q']['min_bid'] = -15.0
    gen, gencost = off2case(gen0, gencost0, offers, bids, lim)

    gen1 = gen0.copy()
    gen1[:, [GEN_STATUS, PMIN, PMAX, QMIN, QMAX]] = array([
        [1,  10, 10, -15,  10],
        [1,  12, 50, -10,  10],
        [1, -10,  0,  -5,   0],
        [1,  12, 25,   0,  30],
        [0, -30,  0,   0, 7.5]
    ])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0[:, :12].copy()
    gencost1[:, NCOST - 1:NCOST + 9] = array([
        [2,   0,     0,  10,   100,   0,    0,  0,    0],
        [3,   0,     0,  20,   500,  50, 2450,  0,    0],
        [3, -30, -2600, -20, -2000,   0,    0,  0,    0],
        [2,   0,     0,  25,  1250,   0,    0,  0,    0],
        [4, -30,     0, -20,  1000, -10, 2000,  0, 3000],
        [4, -15,   150,   0,     0,   5,   50, 10,  150],
        [3, -10,     0,   0,     0,  10,   50,  0,    0],
        [2, -15,   -75,   0,     0,   0,    0,  0,    0],
        [3,   0,     0,  15,    15,  30,  165,  0,    0],
        [2,   0,     0,   0,     0,   0,    0,  0,    0]
    ])
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'PQ offers & PQ bids, for gen, no P, no shutdown';
    gen2 = gen0.copy()
    gen2[0, PMIN] = 0
    offers['P']['qty'] = array([[0, 40], [20, 30], [25, 25]], float)
    gen, gencost = off2case(gen2, gencost0, offers, bids, lim)

    gen1[0, [PMIN, PMAX, QMIN, QMAX]] = array([0, 0, -15, 10])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1[0, NCOST:NCOST + 9] = gencost0[0, NCOST:NCOST + 9]
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'PQ offers & PQ bids, for gen, no Q, no shutdown';
    offers['P']['qty'] = array([[10, 40], [20, 30], [25, 25]], float)
    offers['Q']['qty'] = array([[ 5,  5], [ 0, 10], [15, 15]], float)
    bids['Q']['qty'] = array([15, 0, 15, 15, 0], float)
    gen, gencost = off2case(gen0, gencost0, offers, bids, lim)

    gen1[0, [PMIN, PMAX, QMIN, QMAX]] = array([10, 10, -15, 10])    ## restore original
    gen1[1, [PMIN, PMAX, QMIN, QMAX]] = array([12, 50,   0,  0])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1[ix_([0, 1, 6], range(NCOST, NCOST + 9))] = array([
        [2, 0, 0, 10, 100,  0,    0, 0, 0],
        [3, 0, 0, 20, 500, 50, 2450, 0, 0],
        [2, 0, 0,  0,   0,  0,    0, 0, 0]
    ])
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'PQ offers & PQ bids, lim.P/Q.max_offer/min_bid, multi-block';
    offers['P']['qty'] = array([[10,  40], [20, 30], [25, 25]], float)
    offers['P']['prc'] = array([[10, 100], [25, 65], [50, 90]], float)
    bids['P']['qty'] = array([[10,   0], [12, 18]], float)
    bids['P']['prc'] = array([[100, 60], [70, 10]], float)
    offers['Q']['qty'] = array([[5, 5], [10, 10], [15, 15]], float)
    offers['Q']['prc'] = array([[10, 20], [5, 60], [1, 10]], float)
    bids['Q']['qty'] = array([15, 10, 10, 15, 0], float)
    bids['Q']['prc'] = array([-10, 0, 5, -20, 10], float)
    lim['Q']['max_offer'] = 50.0
    lim['Q']['min_bid'] = -15.0
    gen, gencost = off2case(gen0, gencost0, offers, bids, lim)

    gen1 = gen0.copy()
    gen1[:, [GEN_STATUS, PMIN, PMAX, QMIN, QMAX]] = array([
        [1,  10, 10, -15, 10],
        [1,  12, 50, -10, 10],
        [1, -10,  0,  -5,  0],
        [1,  12, 25,   0, 30],
        [0, -30,  0,   0,  7.5]
    ])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0[:, :12].copy()
    gencost1[:, NCOST:NCOST + 9] = array([
        [2,   0,     0,  10,  100,   0,    0,  0,    0],
        [3,   0,     0,  20,  500,  50, 2450,  0,    0],
        [2, -10, -1000,   0,    0,   0,    0,  0,    0],
        [2,   0,     0,  25, 1250,   0,    0,  0,    0],
        [4, -30,     0, -20, 1000, -10, 2000,  0, 3000],
        [4, -15,   150,   0,    0,   5,   50, 10,  150],
        [3, -10,     0,   0,    0,  10,   50,  0,    0],
        [2, -10,   -50,   0,    0,   0,    0,  0,    0],
        [3,   0,     0,  15,   15,  30,  165,  0,    0],
        [2,   0,     0,   0,    0,   0,    0,  0,    0]
    ])
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'PQ offers & PQ bids, zero Q load w/P bid, shutdown bugfix';
    gen1 = gen0.copy()
    gen1[4, [QG, QMIN, QMAX]] = 0
    gen, gencost = off2case(gen1, gencost0, offers, bids, lim)

    gen1[:, [PMIN, PMAX, QMIN, QMAX]] = array([
        [ 10, 10, -15, 10],
        [ 12, 50, -10, 10],
        [-10,  0,  -5,  0],
        [ 12, 25,   0, 30],
        [-12,  0,   0,  0]
    ])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0[:, :12].copy()
    gencost1[:, NCOST - 1:NCOST + 9] = array([
        [2,   0,     0, 10,  100,  0,    0,  0,   0],
        [3,   0,     0, 20,  500, 50, 2450,  0,   0],
        [2, -10, -1000,  0,    0,  0,    0,  0,   0],
        [2,   0,     0, 25, 1250,  0,    0,  0,   0],
        [2, -12,  -840,  0,    0,  0,    0,  0,   0],
        [4, -15,   150,  0,    0,  5,   50, 10, 150],
        [3, -10,     0,  0,    0, 10,   50,  0,   0],
        [2, -10,   -50,  0,    0,  0,    0,  0,   0],
        [3,   0,     0, 15,   15, 30,  165,  0,   0],
        [2,   0,     0,  0,    0,  0,    0,  0,   0]
    ])
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t = 'PQ offers & PQ bids, non-zero Q load w/no P bid, shutdown bugfix';
    offers['P']['qty'] = array([[10,  40], [20, 30], [25, 25]], float)
    offers['P']['prc'] = array([[10, 100], [25, 65], [50, 90]], float)
    bids['P']['qty'] = array([[0, 10], [12, 18]], float)
    bids['P']['prc'] = array([[100, 40], [70, 10]], float)
    offers['Q']['qty'] = array([[ 5,  5], [10, 10], [15, 15]], float)
    offers['Q']['prc'] = array([[10, 20], [ 5, 60], [ 1, 10]], float)
    bids['Q']['qty'] = array([ 15, 10, 15,  15,  0], float)
    bids['Q']['prc'] = array([-10,  0,  5, -20, 10], float)
    lim['Q']['max_offer'] = 50.0
    lim['Q']['min_bid'] = -15.0
    gen, gencost = off2case(gen0, gencost0, offers, bids, lim)

    gen1 = gen0.copy()
    gen1[:, [GEN_STATUS, PMIN, PMAX, QMIN, QMAX]] = array([
        [1,  10, 10, -15, 10],
        [1,  12, 50, -10, 10],
        [0, -30,  0, -15,  0],
        [1,  12, 25,   0, 30],
        [0, -30,  0,   0, 7.5]
    ])
    t_is( gen, gen1, 8, [t, ' - gen'] )

    gencost1 = gencost0[:, :12].copy()
    gencost1[:, NCOST - 1:NCOST + 9] = array([
        [2,   0,   0,  10,  100,   0,    0,  0,    0],
        [3,   0,   0,  20,  500,  50, 2450,  0,    0],
        [4, -30,   0, -20, 1000, -10, 2000,  0, 3000],
        [2,   0,   0,  25, 1250,   0,    0,  0,    0],
        [4, -30,   0, -20, 1000, -10, 2000,  0, 3000],
        [4, -15, 150,   0,    0,   5,   50, 10,  150],
        [3, -10,   0,   0,    0,  10,   50,  0,    0],
        [3, -20, -15, -10,  -10,   0,    0,  0,    0],
        [3,   0,   0,  15,   15,  30,  165,  0,    0],
        [2,   0,   0,   0,    0,   0,    0,  0,    0]
    ])
    t_is( gencost, gencost1, 8, [t, ' - gencost'] )

    t_end()

Example 21

Project: PYPOWER Source File: t_opf_ipopt.py
def t_opf_ipopt(quiet=False):
    """Tests for IPOPT-based AC optimal power flow.

    @author: Ray Zimmerman (PSERC Cornell)
    """
    num_tests = 101

    t_begin(num_tests, quiet)

    tdir = dirname(__file__)
    casefile = join(tdir, 't_case9_opf')
    verbose = 0#not quiet

    t0 = 'IPOPT : '
    ppopt = ppoption(OPF_VIOLATION=1e-6, PDIPM_GRADTOL=1e-8,
                   PDIPM_COMPTOL=1e-8, PDIPM_COSTTOL=1e-9)
    ppopt = ppoption(ppopt, OUT_ALL=0, VERBOSE=verbose, OPF_ALG=580)

    ## set up indices
    ib_data     = r_[arange(BUS_AREA + 1), arange(BASE_KV, VMIN + 1)]
    ib_voltage  = arange(VM, VA + 1)
    ib_lam      = arange(LAM_P, LAM_Q + 1)
    ib_mu       = arange(MU_VMAX, MU_VMIN + 1)
    ig_data     = r_[[GEN_BUS, QMAX, QMIN], arange(MBASE, APF + 1)]
    ig_disp     = array([PG, QG, VG])
    ig_mu       = arange(MU_PMAX, MU_QMIN + 1)
    ibr_data    = arange(ANGMAX + 1)
    ibr_flow    = arange(PF, QT + 1)
    ibr_mu      = array([MU_SF, MU_ST])
    ibr_angmu   = array([MU_ANGMIN, MU_ANGMAX])

    ## get solved AC power flow case from MAT-file
    soln9_opf = loadmat(join(tdir, 'soln9_opf.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf['bus_soln']
    gen_soln = soln9_opf['gen_soln']
    branch_soln = soln9_opf['branch_soln']
    f_soln = soln9_opf['f_soln'][0]

    ## run OPF
    t = t0
    r = runopf(casefile, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ## run with automatic conversion of single-block pwl to linear costs
    t = ''.join([t0, '(single-block PWL) : '])
    ppc = loadcase(casefile)
    ppc['gencost'][2, NCOST] = 2
    r = runopf(ppc, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    xr = r_[r['var']['val']['Va'], r['var']['val']['Vm'], r['var']['val']['Pg'],
            r['var']['val']['Qg'], 0, r['var']['val']['y']]
    t_is(r['x'], xr, 8, [t, 'check on raw x returned from OPF'])

    ## get solved AC power flow case from MAT-file
    soln9_opf_Plim = loadmat(join(tdir, 'soln9_opf_Plim.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_Plim['bus_soln']
    gen_soln = soln9_opf_Plim['gen_soln']
    branch_soln = soln9_opf_Plim['branch_soln']
    f_soln = soln9_opf_Plim['f_soln'][0]

    ## run OPF with active power line limits
    t = ''.join([t0, '(P line lim) : '])
    ppopt1 = ppoption(ppopt, OPF_FLOW_LIM=1)
    r = runopf(casefile, ppopt1)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ##-----  test OPF with quadratic gen costs moved to generalized costs  -----
    ppc = loadcase(casefile)
    ppc['gencost'] = array([
        [2,   1500, 0,   3,   0.11,    5,   0],
        [2,   2000, 0,   3,   0.085,   1.2, 0],
        [2,   3000, 0,   3,   0.1225,  1,   0]
    ])
    r = runopf(ppc, ppopt)
    bus_soln, gen_soln, branch_soln, f_soln, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    branch_soln = branch_soln[:, :MU_ST + 1]

    A = None
    l = array([])
    u = array([])
    nb = ppc['bus'].shape[0]      # number of buses
    ng = ppc['gen'].shape[0]      # number of gens
    thbas = 0;                thend    = thbas + nb
    vbas     = thend;     vend     = vbas + nb
    pgbas    = vend;      pgend    = pgbas + ng
#    qgbas    = pgend;     qgend    = qgbas + ng
    nxyz = 2 * nb + 2 * ng
    N = sparse((ppc['baseMVA'] * ones(ng), (arange(ng), arange(pgbas, pgend))), (ng, nxyz))
    fparm = ones((ng, 1)) * array([[1, 0, 0, 1]])
    ix = argsort(ppc['gen'][:, 0])
    H = 2 * spdiags(ppc['gencost'][ix, 4], 0, ng, ng, 'csr')
    Cw = ppc['gencost'][ix, 5]
    ppc['gencost'][:, 4:7] = 0

    ## run OPF with quadratic gen costs moved to generalized costs
    t = ''.join([t0, 'w/quadratic generalized gen cost : '])
    r = opf(ppc, A, l, u, ppopt, N, fparm, H, Cw)
    f, bus, gen, branch, success = \
            r['f'], r['bus'], r['gen'], r['branch'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(r['cost']['usr'], f, 12, [t, 'user cost'])

    ##-----  run OPF with extra linear user constraints & costs  -----
    ## single new z variable constrained to be greater than or equal to
    ## deviation from 1 pu voltage at bus 1, linear cost on this z
    ## get solved AC power flow case from MAT-file
    soln9_opf_extras1 = loadmat(join(tdir, 'soln9_opf_extras1.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_extras1['bus_soln']
    gen_soln = soln9_opf_extras1['gen_soln']
    branch_soln = soln9_opf_extras1['branch_soln']
    f_soln = soln9_opf_extras1['f_soln'][0]

    row = [0, 0, 1, 1]
    col = [9, 24, 9, 24]
    A = sparse(([-1, 1, 1, 1], (row, col)), (2, 25))
    u = array([Inf, Inf])
    l = array([-1, 1])

    N = sparse(([1], ([0], [24])), (1, 25))    ## new z variable only
    fparm = array([[1, 0, 0, 1]])              ## w = r = z
    H = sparse((1, 1))                ## no quadratic term
    Cw = array([100.0])

    t = ''.join([t0, 'w/extra constraints & costs 1 : '])
    r = opf(casefile, A, l, u, ppopt, N, fparm, H, Cw)
    f, bus, gen, branch, success = \
            r['f'], r['bus'], r['gen'], r['branch'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(r['var']['val']['z'], 0.025419, 6, [t, 'user variable'])
    t_is(r['cost']['usr'], 2.5419, 4, [t, 'user cost'])

    ##-----  test OPF with capability curves  -----
    ppc = loadcase(join(tdir, 't_case9_opfv2'))
    ## remove angle diff limits
    ppc['branch'][0, ANGMAX] =  360
    ppc['branch'][8, ANGMIN] = -360

    ## get solved AC power flow case from MAT-file
    soln9_opf_PQcap = loadmat(join(tdir, 'soln9_opf_PQcap.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_PQcap['bus_soln']
    gen_soln = soln9_opf_PQcap['gen_soln']
    branch_soln = soln9_opf_PQcap['branch_soln']
    f_soln = soln9_opf_PQcap['f_soln'][0]

    ## run OPF with capability curves
    t = ''.join([t0, 'w/capability curves : '])
    r = runopf(ppc, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ##-----  test OPF with angle difference limits  -----
    ppc = loadcase(join(tdir, 't_case9_opfv2'))
    ## remove capability curves
    ppc['gen'][ix_(arange(1, 3),
                   [PC1, PC2, QC1MIN, QC1MAX, QC2MIN, QC2MAX])] = zeros((2, 6))

    ## get solved AC power flow case from MAT-file
    soln9_opf_ang = loadmat(join(tdir, 'soln9_opf_ang.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_ang['bus_soln']
    gen_soln = soln9_opf_ang['gen_soln']
    branch_soln = soln9_opf_ang['branch_soln']
    f_soln = soln9_opf_ang['f_soln'][0]

    ## run OPF with angle difference limits
    t = ''.join([t0, 'w/angle difference limits : '])
    r = runopf(ppc, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  1, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(branch[:, ibr_angmu ], branch_soln[:, ibr_angmu ],  2, [t, 'branch angle mu'])

    ##-----  test OPF with ignored angle difference limits  -----
    ## get solved AC power flow case from MAT-file
    soln9_opf = loadmat(join(tdir, 'soln9_opf.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf['bus_soln']
    gen_soln = soln9_opf['gen_soln']
    branch_soln = soln9_opf['branch_soln']
    f_soln = soln9_opf['f_soln'][0]

    ## run OPF with ignored angle difference limits
    t = ''.join([t0, 'w/ignored angle difference limits : '])
    ppopt1 = ppoption(ppopt, OPF_IGNORE_ANG_LIM=1)
    r = runopf(ppc, ppopt1)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    ## ang limits are not in this solution data, so let's remove them
    branch[0, ANGMAX] =  360
    branch[8, ANGMIN] = -360
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    t_end()

Example 22

Project: PYPOWER Source File: t_opf_pips.py
def t_opf_pips(quiet=False):
    """Tests for PIPS-based AC optimal power flow.

    @author: Ray Zimmerman (PSERC Cornell)
    """
    num_tests = 101

    t_begin(num_tests, quiet)

    tdir = dirname(__file__)
    casefile = join(tdir, 't_case9_opf')
    verbose = 0#not quiet

    t0 = 'PIPS : '
    ppopt = ppoption(OPF_VIOLATION=1e-6, PDIPM_GRADTOL=1e-8,
                   PDIPM_COMPTOL=1e-8, PDIPM_COSTTOL=1e-9)
    ppopt = ppoption(ppopt, OUT_ALL=0, VERBOSE=verbose, OPF_ALG=560)

    ## set up indices
    ib_data     = r_[arange(BUS_AREA + 1), arange(BASE_KV, VMIN + 1)]
    ib_voltage  = arange(VM, VA + 1)
    ib_lam      = arange(LAM_P, LAM_Q + 1)
    ib_mu       = arange(MU_VMAX, MU_VMIN + 1)
    ig_data     = r_[[GEN_BUS, QMAX, QMIN], arange(MBASE, APF + 1)]
    ig_disp     = array([PG, QG, VG])
    ig_mu       = arange(MU_PMAX, MU_QMIN + 1)
    ibr_data    = arange(ANGMAX + 1)
    ibr_flow    = arange(PF, QT + 1)
    ibr_mu      = array([MU_SF, MU_ST])
    ibr_angmu   = array([MU_ANGMIN, MU_ANGMAX])

    ## get solved AC power flow case from MAT-file
    soln9_opf = loadmat(join(tdir, 'soln9_opf.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf['bus_soln']
    gen_soln = soln9_opf['gen_soln']
    branch_soln = soln9_opf['branch_soln']
    f_soln = soln9_opf['f_soln'][0]

    ## run OPF
    t = t0
    r = runopf(casefile, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ## run with automatic conversion of single-block pwl to linear costs
    t = ''.join([t0, '(single-block PWL) : '])
    ppc = loadcase(casefile)
    ppc['gencost'][2, NCOST] = 2
    r = runopf(ppc, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    xr = r_[r['var']['val']['Va'], r['var']['val']['Vm'], r['var']['val']['Pg'],
            r['var']['val']['Qg'], 0, r['var']['val']['y']]
    t_is(r['x'], xr, 8, [t, 'check on raw x returned from OPF'])

    ## get solved AC power flow case from MAT-file
    soln9_opf_Plim = loadmat(join(tdir, 'soln9_opf_Plim.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_Plim['bus_soln']
    gen_soln = soln9_opf_Plim['gen_soln']
    branch_soln = soln9_opf_Plim['branch_soln']
    f_soln = soln9_opf_Plim['f_soln'][0]

    ## run OPF with active power line limits
    t = ''.join([t0, '(P line lim) : '])
    ppopt1 = ppoption(ppopt, OPF_FLOW_LIM=1)
    r = runopf(casefile, ppopt1)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ##-----  test OPF with quadratic gen costs moved to generalized costs  -----
    ppc = loadcase(casefile)
    ppc['gencost'] = array([
        [2,   1500, 0,   3,   0.11,    5,   0],
        [2,   2000, 0,   3,   0.085,   1.2, 0],
        [2,   3000, 0,   3,   0.1225,  1,   0]
    ])
    r = runopf(ppc, ppopt)
    bus_soln, gen_soln, branch_soln, f_soln, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    branch_soln = branch_soln[:, :MU_ST + 1]

    A = None
    l = array([])
    u = array([])
    nb = ppc['bus'].shape[0]      # number of buses
    ng = ppc['gen'].shape[0]      # number of gens
    thbas = 0;                thend    = thbas + nb
    vbas     = thend;     vend     = vbas + nb
    pgbas    = vend;      pgend    = pgbas + ng
#    qgbas    = pgend;     qgend    = qgbas + ng
    nxyz = 2 * nb + 2 * ng
    N = sparse((ppc['baseMVA'] * ones(ng), (arange(ng), arange(pgbas, pgend))), (ng, nxyz))
    fparm = ones((ng, 1)) * array([[1, 0, 0, 1]])
    ix = argsort(ppc['gen'][:, 0])
    H = 2 * spdiags(ppc['gencost'][ix, 4], 0, ng, ng, 'csr')
    Cw = ppc['gencost'][ix, 5]
    ppc['gencost'][:, 4:7] = 0

    ## run OPF with quadratic gen costs moved to generalized costs
    t = ''.join([t0, 'w/quadratic generalized gen cost : '])
    r = opf(ppc, A, l, u, ppopt, N, fparm, H, Cw)
    f, bus, gen, branch, success = \
            r['f'], r['bus'], r['gen'], r['branch'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(r['cost']['usr'], f, 12, [t, 'user cost'])

    ##-----  run OPF with extra linear user constraints & costs  -----
    ## single new z variable constrained to be greater than or equal to
    ## deviation from 1 pu voltage at bus 1, linear cost on this z
    ## get solved AC power flow case from MAT-file
    soln9_opf_extras1 = loadmat(join(tdir, 'soln9_opf_extras1.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_extras1['bus_soln']
    gen_soln = soln9_opf_extras1['gen_soln']
    branch_soln = soln9_opf_extras1['branch_soln']
    f_soln = soln9_opf_extras1['f_soln'][0]

    row = [0, 0, 1, 1]
    col = [9, 24, 9, 24]
    A = sparse(([-1, 1, 1, 1], (row, col)), (2, 25))
    u = array([Inf, Inf])
    l = array([-1, 1])

    N = sparse(([1], ([0], [24])), (1, 25))    ## new z variable only
    fparm = array([[1, 0, 0, 1]])              ## w = r = z
    H = sparse((1, 1))                ## no quadratic term
    Cw = array([100.0])

    t = ''.join([t0, 'w/extra constraints & costs 1 : '])
    r = opf(casefile, A, l, u, ppopt, N, fparm, H, Cw)
    f, bus, gen, branch, success = \
            r['f'], r['bus'], r['gen'], r['branch'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(r['var']['val']['z'], 0.025419, 6, [t, 'user variable'])
    t_is(r['cost']['usr'], 2.5419, 4, [t, 'user cost'])

    ##-----  test OPF with capability curves  -----
    ppc = loadcase(join(tdir, 't_case9_opfv2'))
    ## remove angle diff limits
    ppc['branch'][0, ANGMAX] =  360
    ppc['branch'][8, ANGMIN] = -360

    ## get solved AC power flow case from MAT-file
    soln9_opf_PQcap = loadmat(join(tdir, 'soln9_opf_PQcap.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_PQcap['bus_soln']
    gen_soln = soln9_opf_PQcap['gen_soln']
    branch_soln = soln9_opf_PQcap['branch_soln']
    f_soln = soln9_opf_PQcap['f_soln'][0]

    ## run OPF with capability curves
    t = ''.join([t0, 'w/capability curves : '])
    r = runopf(ppc, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ##-----  test OPF with angle difference limits  -----
    ppc = loadcase(join(tdir, 't_case9_opfv2'))
    ## remove capability curves
    ppc['gen'][ix_(arange(1, 3),
                   [PC1, PC2, QC1MIN, QC1MAX, QC2MIN, QC2MAX])] = zeros((2, 6))

    ## get solved AC power flow case from MAT-file
    soln9_opf_ang = loadmat(join(tdir, 'soln9_opf_ang.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_ang['bus_soln']
    gen_soln = soln9_opf_ang['gen_soln']
    branch_soln = soln9_opf_ang['branch_soln']
    f_soln = soln9_opf_ang['f_soln'][0]

    ## run OPF with angle difference limits
    t = ''.join([t0, 'w/angle difference limits : '])
    r = runopf(ppc, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  1, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(branch[:, ibr_angmu ], branch_soln[:, ibr_angmu ],  2, [t, 'branch angle mu'])

    ##-----  test OPF with ignored angle difference limits  -----
    ## get solved AC power flow case from MAT-file
    soln9_opf = loadmat(join(tdir, 'soln9_opf.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf['bus_soln']
    gen_soln = soln9_opf['gen_soln']
    branch_soln = soln9_opf['branch_soln']
    f_soln = soln9_opf['f_soln'][0]

    ## run OPF with ignored angle difference limits
    t = ''.join([t0, 'w/ignored angle difference limits : '])
    ppopt1 = ppoption(ppopt, OPF_IGNORE_ANG_LIM=1)
    r = runopf(ppc, ppopt1)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    ## ang limits are not in this solution data, so let's remove them
    branch[0, ANGMAX] =  360
    branch[8, ANGMIN] = -360
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    t_end()

Example 23

Project: PYPOWER Source File: t_opf_pips_sc.py
def t_opf_pips_sc(quiet=False):
    """Tests for step-controlled PIPS-based AC optimal power flow.

    @author: Ray Zimmerman (PSERC Cornell)
    """
    num_tests = 101

    t_begin(num_tests, quiet)

    tdir = dirname(__file__)
    casefile = join(tdir, 't_case9_opf')
    verbose = 0#not quiet

    t0 = 'PIPS-sc : '
    ppopt = ppoption(OPF_VIOLATION=1e-6, PDIPM_GRADTOL=1e-8,
                     PDIPM_COMPTOL=1e-8, PDIPM_COSTTOL=1e-9)
    ppopt = ppoption(ppopt, OUT_ALL=0, VERBOSE=verbose, OPF_ALG=565)

    ## set up indices
    ib_data     = r_[arange(BUS_AREA + 1), arange(BASE_KV, VMIN + 1)]
    ib_voltage  = arange(VM, VA + 1)
    ib_lam      = arange(LAM_P, LAM_Q + 1)
    ib_mu       = arange(MU_VMAX, MU_VMIN + 1)
    ig_data     = r_[[GEN_BUS, QMAX, QMIN], arange(MBASE, APF + 1)]
    ig_disp     = array([PG, QG, VG])
    ig_mu       = arange(MU_PMAX, MU_QMIN + 1)
    ibr_data    = arange(ANGMAX + 1)
    ibr_flow    = arange(PF, QT + 1)
    ibr_mu      = array([MU_SF, MU_ST])
    ibr_angmu   = array([MU_ANGMIN, MU_ANGMAX])

    ## get solved AC power flow case from MAT-file
    soln9_opf = loadmat(join(tdir, 'soln9_opf.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf['bus_soln']
    gen_soln = soln9_opf['gen_soln']
    branch_soln = soln9_opf['branch_soln']
    f_soln = soln9_opf['f_soln'][0]

    ## run OPF
    t = t0
    r = runopf(casefile, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ## run with automatic conversion of single-block pwl to linear costs
    t = ''.join([t0, '(single-block PWL) : '])
    ppc = loadcase(casefile)
    ppc['gencost'][2, NCOST] = 2
    r = runopf(ppc, ppopt)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    xr = r_[r['var']['val']['Va'], r['var']['val']['Vm'], r['var']['val']['Pg'],
            r['var']['val']['Qg'], 0, r['var']['val']['y']]
    t_is(r['x'], xr, 8, [t, 'check on raw x returned from OPF'])

    ## get solved AC power flow case from MAT-file
    soln9_opf_Plim = loadmat(join(tdir, 'soln9_opf_Plim.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_Plim['bus_soln']
    gen_soln = soln9_opf_Plim['gen_soln']
    branch_soln = soln9_opf_Plim['branch_soln']
    f_soln = soln9_opf_Plim['f_soln'][0]

    ## run OPF with active power line limits
    t = ''.join([t0, '(P line lim) : '])
    ppopt1 = ppoption(ppopt, OPF_FLOW_LIM=1)
    r = runopf(casefile, ppopt1)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ##-----  test OPF with quadratic gen costs moved to generalized costs  -----
    ppc = loadcase(casefile)
    ppc['gencost'] = array([
        [2,   1500, 0,   3,   0.11,    5,   0],
        [2,   2000, 0,   3,   0.085,   1.2, 0],
        [2,   3000, 0,   3,   0.1225,  1,   0]
    ])
    r = runopf(ppc, ppopt)
    bus_soln, gen_soln, branch_soln, f_soln, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    branch_soln = branch_soln[:, :MU_ST + 1]

    A = None
    l = array([])
    u = array([])
    nb = ppc['bus'].shape[0]      # number of buses
    ng = ppc['gen'].shape[0]      # number of gens
    thbas = 0;            thend    = thbas + nb
    vbas     = thend;     vend     = vbas + nb
    pgbas    = vend;      pgend    = pgbas + ng
#    qgbas    = pgend;     qgend    = qgbas + ng
    nxyz = 2 * nb + 2 * ng
    N = sparse((ppc['baseMVA'] * ones(ng), (arange(ng), arange(pgbas, pgend))), (ng, nxyz))
    fparm = ones((ng, 1)) * array([[1, 0, 0, 1]])
    ix = argsort(ppc['gen'][:, 0])
    H = 2 * spdiags(ppc['gencost'][ix, 4], 0, ng, ng, 'csr')
    Cw = ppc['gencost'][ix, 5]
    ppc['gencost'][:, 4:7] = 0

    ## run OPF with quadratic gen costs moved to generalized costs
    t = ''.join([t0, 'w/quadratic generalized gen cost : '])
    r = opf(ppc, A, l, u, ppopt, N, fparm, H, Cw)
    f, bus, gen, branch, success = \
            r['f'], r['bus'], r['gen'], r['branch'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(r['cost']['usr'], f, 12, [t, 'user cost'])

    ##-----  run OPF with extra linear user constraints & costs  -----
    ## single new z variable constrained to be greater than or equal to
    ## deviation from 1 pu voltage at bus 1, linear cost on this z
    ## get solved AC power flow case from MAT-file
    soln9_opf_extras1 = loadmat(join(tdir, 'soln9_opf_extras1.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_extras1['bus_soln']
    gen_soln = soln9_opf_extras1['gen_soln']
    branch_soln = soln9_opf_extras1['branch_soln']
    f_soln = soln9_opf_extras1['f_soln'][0]

    row = [0, 0, 1, 1]
    col = [9, 24, 9, 24]
    A = sparse(([-1, 1, 1, 1], (row, col)), (2, 25))
    u = array([Inf, Inf])
    l = array([-1, 1])

    N = sparse(([1], ([0], [24])), (1, 25))    ## new z variable only
    fparm = array([[1, 0, 0, 1]])              ## w = r = z
    H = sparse((1, 1))                ## no quadratic term
    Cw = array([100.0])

    t = ''.join([t0, 'w/extra constraints & costs 1 : '])
    r = opf(casefile, A, l, u, ppopt, N, fparm, H, Cw)
    f, bus, gen, branch, success = \
            r['f'], r['bus'], r['gen'], r['branch'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(r['var']['val']['z'], 0.025419, 6, [t, 'user variable'])
    t_is(r['cost']['usr'], 2.5419, 4, [t, 'user cost'])

    ##-----  test OPF with capability curves  -----
    ppc = loadcase(join(tdir, 't_case9_opfv2'))
    ## remove angle diff limits
    ppc['branch'][0, ANGMAX] = 360
    ppc['branch'][8, ANGMIN] = -360

    ## get solved AC power flow case from MAT-file
    soln9_opf_PQcap = loadmat(join(tdir, 'soln9_opf_PQcap.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_PQcap['bus_soln']
    gen_soln = soln9_opf_PQcap['gen_soln']
    branch_soln = soln9_opf_PQcap['branch_soln']
    f_soln = soln9_opf_PQcap['f_soln'][0]

    ## run OPF with capability curves
    t = ''.join([t0, 'w/capability curves : '])
    r = runopf(ppc, ppopt)
    f, bus, gen, branch, success = \
            r['f'], r['bus'], r['gen'], r['branch'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    ##-----  test OPF with angle difference limits  -----
    ppc = loadcase(join(tdir, 't_case9_opfv2'))
    ## remove capability curves
    ppc['gen'][ix_(arange(1, 3),
                   [PC1, PC2, QC1MIN, QC1MAX, QC2MIN, QC2MAX])] = zeros((2, 6))

    ## get solved AC power flow case from MAT-file
    soln9_opf_ang = loadmat(join(tdir, 'soln9_opf_ang.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf_ang['bus_soln']
    gen_soln = soln9_opf_ang['gen_soln']
    branch_soln = soln9_opf_ang['branch_soln']
    f_soln = soln9_opf_ang['f_soln'][0]

    ## run OPF with angle difference limits
    t = ''.join([t0, 'w/angle difference limits : '])
    r = runopf(ppc, ppopt)
    f, bus, gen, branch, success = \
            r['f'], r['bus'], r['gen'], r['branch'], r['success']
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  1, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])
    t_is(branch[:, ibr_angmu ], branch_soln[:, ibr_angmu ],  2, [t, 'branch angle mu'])

    ##-----  test OPF with ignored angle difference limits  -----
    ## get solved AC power flow case from MAT-file
    soln9_opf = loadmat(join(tdir, 'soln9_opf.mat'), struct_as_record=True)
    ## defines bus_soln, gen_soln, branch_soln, f_soln
    bus_soln = soln9_opf['bus_soln']
    gen_soln = soln9_opf['gen_soln']
    branch_soln = soln9_opf['branch_soln']
    f_soln = soln9_opf['f_soln'][0]

    ## run OPF with ignored angle difference limits
    t = ''.join([t0, 'w/ignored angle difference limits : '])
    ppopt1 = ppoption(ppopt, OPF_IGNORE_ANG_LIM=1)
    r = runopf(ppc, ppopt1)
    bus, gen, branch, f, success = \
            r['bus'], r['gen'], r['branch'], r['f'], r['success']
    ## ang limits are not in this solution data, so let's remove them
    branch[0, ANGMAX] =  360
    branch[8, ANGMIN] = -360
    t_ok(success, [t, 'success'])
    t_is(f, f_soln, 3, [t, 'f'])
    t_is(   bus[:, ib_data   ],    bus_soln[:, ib_data   ], 10, [t, 'bus data'])
    t_is(   bus[:, ib_voltage],    bus_soln[:, ib_voltage],  3, [t, 'bus voltage'])
    t_is(   bus[:, ib_lam    ],    bus_soln[:, ib_lam    ],  3, [t, 'bus lambda'])
    t_is(   bus[:, ib_mu     ],    bus_soln[:, ib_mu     ],  2, [t, 'bus mu'])
    t_is(   gen[:, ig_data   ],    gen_soln[:, ig_data   ], 10, [t, 'gen data'])
    t_is(   gen[:, ig_disp   ],    gen_soln[:, ig_disp   ],  3, [t, 'gen dispatch'])
    t_is(   gen[:, ig_mu     ],    gen_soln[:, ig_mu     ],  3, [t, 'gen mu'])
    t_is(branch[:, ibr_data  ], branch_soln[:, ibr_data  ], 10, [t, 'branch data'])
    t_is(branch[:, ibr_flow  ], branch_soln[:, ibr_flow  ],  3, [t, 'branch flow'])
    t_is(branch[:, ibr_mu    ], branch_soln[:, ibr_mu    ],  2, [t, 'branch mu'])

    t_end()

Example 24

Project: scikit-learn Source File: metaestimators.py
Function: safe_split
def _safe_split(estimator, X, y, indices, train_indices=None):
    """Create subset of dataset and properly handle kernels."""
    from ..gaussian_process.kernels import Kernel as GPKernel

    if (hasattr(estimator, 'kernel') and callable(estimator.kernel) and
            not isinstance(estimator.kernel, GPKernel)):
        # cannot compute the kernel values with custom function
        raise ValueError("Cannot use a custom kernel function. "
                         "Precompute the kernel matrix instead.")

    if not hasattr(X, "shape"):
        if getattr(estimator, "_pairwise", False):
            raise ValueError("Precomputed kernels or affinity matrices have "
                             "to be passed as arrays or sparse matrices.")
        X_subset = [X[index] for index in indices]
    else:
        if getattr(estimator, "_pairwise", False):
            # X is a precomputed square kernel matrix
            if X.shape[0] != X.shape[1]:
                raise ValueError("X should be a square kernel matrix")
            if train_indices is None:
                X_subset = X[np.ix_(indices, indices)]
            else:
                X_subset = X[np.ix_(indices, train_indices)]
        else:
            X_subset = safe_indexing(X, indices)

    if y is not None:
        y_subset = safe_indexing(y, indices)
    else:
        y_subset = None

    return X_subset, y_subset

Example 25

Project: pebl Source File: data.py
Function: sub_set
    def subset(self, variables=None, samples=None):
        """Returns a subset of the dataset (and metadata).
        
        Specify the variables and samples for creating a subset of the data.
        variables and samples should be a list of ids. If not specified, it is
        assumed to be all variables or samples. 

        Some examples:
        
            - d.subset([3], [4])
            - d.subset([3,1,2])
            - d.subset(samples=[5,2,7,1])
        
        Note: order matters! d.subset([3,1,2]) != d.subset([1,2,3])

        """

        variables = variables if variables is not None else range(self.variables.size)
        samples = samples if samples is not None else range(self.samples.size)
        skip_stats = not (self.has_interventions or self.has_missing)
        d = Dataset(
            self.observations[N.ix_(samples,variables)],
            self.missing[N.ix_(samples,variables)],
            self.interventions[N.ix_(samples,variables)],
            self.variables[variables],
            self.samples[samples],
            skip_stats = skip_stats
        )
        
        # if self does not have interventions or missing, the subset can't.
        if skip_stats:
            d._has_interventions = False
            d._has_missing = False

        return d

Example 26

Project: pyphi Source File: concept.py
    def _relevant_connections(self, subsystem):
        """Identify connections that “matter” to this concept.

        For a core cause, the important connections are those which connect the
        purview to the mechanism; for a core effect they are the connections
        from the mechanism to the purview.

        Returns an |n x n| matrix, where `n` is the number of nodes in this
        corresponding subsystem, that identifies connections that “matter” to
        this MICE:

        ``direction == 'past'``:
            ``relevant_connections[i,j]`` is ``1`` if node ``i`` is in the
            cause purview and node ``j`` is in the mechanism (and ``0``
            otherwise).

        ``direction == 'future'``:
            ``relevant_connections[i,j]`` is ``1`` if node ``i`` is in the
            mechanism and node ``j`` is in the effect purview (and ``0``
            otherwise).

        Args:
            subsystem (Subsystem): The subsystem of this mice

        Returns:
            cm (np.ndarray): A |n x n| matrix of connections, where `n` is the
                size of the subsystem.
        """
        if self.direction == DIRECTIONS[PAST]:
            _from, to = self.purview, self.mechanism
        elif self.direction == DIRECTIONS[FUTURE]:
            _from, to = self.mechanism, self.purview

        cm = utils.relevant_connections(subsystem.network.size, _from, to)
        # Submatrix for this subsystem's nodes
        return cm[np.ix_(subsystem.node_indices, subsystem.node_indices)]

Example 27

Project: pyphi Source File: utils.py
Function: fully_connected
def fully_connected(cm, nodes1, nodes2):
    """Test connectivity of one set of nodes to another.

    Args:
        cm (``np.ndarrray``): The connectivity matrix
        nodes1 (tuple[int]): The nodes whose outputs to ``nodes2`` will be
            tested.
        nodes2 (tuple[int]): The nodes whose inputs from ``nodes1`` will
            be tested.

    Returns:
        bool: Returns True if all elements in ``nodes1`` output to
            some element in ``nodes2`` AND all elements in ``nodes2``
            have an input from some element in ``nodes1``. Otherwise
            return False. Return True if either set of nodes is empty.
    """
    if not nodes1 or not nodes2:
        return True

    cm = cm[np.ix_(nodes1, nodes2)]

    # Do all nodes have at least one connection?
    return cm.sum(0).all() and cm.sum(1).all()

Example 28

Project: aracna Source File: ExternalStrategy.py
    def _getNext(self):
        '''Get the next point to try.  This reads from the file
        self.motionFile'''

        #print 'TRYING READ STDOUT'
        #stdout = self.proc.stdout.read()
        #print 'TRYING READ STDERR'
        #stderr = self.proc.stderr.read()

        #print 'STDOUT:'
        #print stdout
        #print 'STDERR:'
        #print stderr

        if self.orgId == self.generationSize:
            print 'Restarting process after %d runs, push enter when ready...' % self.generationSize
            raw_input()
            #sleep(10)
            if self.spawnProc:
                print 'Continuing with neatfile', self.nextNeatFile
                self.proc = Process((self.executable,
                                     '-O', self.prefix, '-R', '102', '-I',
                                     self.datFile, '-X', self.nextNeatFile, '-XG', '1'))
            self.genId += 1
            self.orgId = 0

        #print 'On iteration', self.orgId+1

        while True:

            if self.spawnProc:
                out = self.proc.read()
                if out != '':
                    #print 'Got stdout:'
                    #print out
                    pass
                out = self.proc.readerr()
                if out != '':
                    #print 'Got stderr:'
                    #print out
                    pass

            try:
                ff = open(self.motionFile, 'r')
            except IOError:
                print 'File does not exist yet'
                sleep(1)
                continue

            lines = ff.readlines()
            nLines = len(lines)
            if nLines < self.expectedLines:
                print '   only %d of %d lines, waiting...' % (nLines, self.expectedLines)
                ff.close()
                sleep(.5)
                continue
            break

        self.orgId += 1

        for ii,line in enumerate(lines[self.junkPoints:]):
            #print 'line', ii, 'is', line
            nums = [float(xx) for xx in line.split()]
            if ii == 0:
                rawPositions = array(nums)
            else:
                rawPositions = vstack((rawPositions, array(nums)))

        # swap and scale rawPositions appropriately
        rawPositions = rawPositions.T[ix_(self.motorColumns)].T  # remap to right columns

        #print 'First few lines of rawPositions are:'
        #for ii in range(10):
        #    print prettyVec(rawPositions[ii,:], prec=2)
        
        # Average over every self.avgPoints
        for ii in range(self.expectedLines/self.avgPoints):
            temp = mean(rawPositions[ii:(ii+self.avgPoints),:], 0)
            if ii == 0:
                positions = temp
            else:
                positions = vstack((positions, temp))

        #print 'First few lines of positions are:'
        #for ii in range(10):
        #    print prettyVec(positions[ii,:], prec=2)
        
        # scale from [-1,1] to [0,1]
        positions += 1
        positions *= .5
        # scale from [0,1] to appropriate ranges
        innerIdx = [0, 2, 4, 6]
        outerIdx = [1, 3, 5, 7]
        centerIdx = [8]
        for ii in innerIdx:
            positions[:,ii] *= (MAX_INNER - MIN_INNER)
            positions[:,ii] += MIN_INNER
        for ii in outerIdx:
            positions[:,ii] *= (MAX_OUTER - MIN_OUTER)
            positions[:,ii] += MIN_OUTER
        for ii in centerIdx:
            positions[:,ii] *= (MAX_CENTER - MIN_CENTER)
            positions[:,ii] += MIN_CENTER
            
        # append a column of 512s for center
        #positions = hstack((positions,
        #                    NORM_CENTER * ones((positions.shape[0],1))))
        times = linspace(0,12,positions.shape[0])

        # Dump both raw positions and positions to file
        thisIdentifier = '%s_%05d_%03d' % (self.identifier, self.genId, self.orgId)

        ff = open('%s_raw' % thisIdentifier, 'w')
        writeArray(ff, rawPositions)
        ff.close()
        ff = open('%s_filt' % thisIdentifier, 'w')
        writeArray(ff, positions)
        ff.close()
        
        # return function of time
        return lambda time: matInterp(time, times, positions), thisIdentifier

Example 29

Project: diogenes Source File: experiment.py
Function: test_m
    def __test_M(self):
        if self.M_test is not None:
            return self.M_test
        return self.M[np.ix_(self.test_indices, self.sub_col_inds)]

Example 30

Project: diogenes Source File: experiment.py
Function: run
    def run(self):
        """Run the Trial"""
        if self.has_run():
            return self.runs
        runs = []
        for subset, sub_col_names, subset_note in self.subset(
                self.labels, 
                self.col_names, 
                **self.subset_params):
            runs_this_subset = []
            labels_sub = self.labels[subset]
            sub_col_inds = self.__indices_of_sublist(
                self.col_names, 
                sub_col_names)
            # np.ix_ from http://stackoverflow.com/questions/30176268/error-when-indexing-with-2-dimensions-in-numpy
            M_sub = self.M[np.ix_(subset, sub_col_inds)]
            cv_cls = self.cv
            cv_kwargs = copy.deepcopy(self.cv_params)
            expected_cv_kwargs = inspect.getargspec(cv_cls.__init__).args
            if 'n' in expected_cv_kwargs:
                cv_kwargs['n'] = labels_sub.shape[0]
            if 'y' in expected_cv_kwargs:
                cv_kwargs['y'] = labels_sub
            if 'labels' in expected_cv_kwargs:
                cv_kwargs['labels'] = labels_sub
            if 'M' in expected_cv_kwargs:
                cv_kwargs['M'] = M_sub
            if 'col_names' in expected_cv_kwargs:
                cv_kwargs['col_names'] = sub_col_names
            cv_inst = cv_cls(**cv_kwargs)
            for fold_idx, (train, test) in enumerate(cv_inst):
                if hasattr(cv_inst, 'cv_note'):
                    cv_note = cv_inst.cv_note()
                else:
                    cv_note = {'fold': fold_idx}
                clf_inst = self.clf(**self.clf_params)
                clf_inst.fit(M_sub[train], labels_sub[train])
                test_indices = subset[test]
                train_indices = subset[train]
                runs_this_subset.append(Run(
                    self.M, 
                    self.labels, 
                    self.col_names, 
                    clf_inst, 
                    train_indices, 
                    test_indices,
                    sub_col_names, 
                    sub_col_inds,
                    subset_note, 
                    cv_note))
            runs.append(runs_this_subset)    
        self.runs = runs
        return self

Example 31

Project: xppy Source File: allinfo.py
    def getFlippedBranch(self, nr, getParts=False):
        '''
        Function returns flipped data of branch number nr;
        additionally it cen return list with indexes of parts
        of the branch (stability wise)
        '''
        (b,p) = self.getBranch(nr,True)
        if b == None:
            return None
            
        # Find starting point
        i = np.ix_(b[:,2]==b[0,2],b[:,3]==b[0,3])[0].ravel()
        # If the branch can't be flipped, return
        try:
            i = i[1]
        except IndexError:
            if getParts:
                return (b,p)
            else:
                return b
        # Change the stability of the first point to the proper one
        b[0,0] = b[1,0]
        # Stack two arrays togather in better order
        retArr = np.vstack((np.flipud(b[0:i,:]),b[i+1:,:]))
        # If we want additional parts array
        if getParts:
            parts = self.findParts(retArr)
            return (retArr, parts)
        # Else, just return the branch array
        return retArr          

Example 32

Project: mne-python Source File: test_cov.py
@slow_test
def test_rank():
    """Test cov rank estimation."""
    # Test that our rank estimation works properly on a simple case
    evoked = read_evokeds(ave_fname, condition=0, baseline=(None, 0),
                          proj=False)
    cov = read_cov(cov_fname)
    ch_names = [ch for ch in evoked.info['ch_names'] if '053' not in ch and
                ch.startswith('EEG')]
    cov = prepare_noise_cov(cov, evoked.info, ch_names, None)
    assert_equal(cov['eig'][0], 0.)  # avg projector should set this to zero
    assert_true((cov['eig'][1:] > 0).all())  # all else should be > 0

    # Now do some more comprehensive tests
    raw_sample = read_raw_fif(raw_fname)

    raw_sss = read_raw_fif(hp_fif_fname)
    raw_sss.add_proj(compute_proj_raw(raw_sss))

    cov_sample = compute_raw_covariance(raw_sample)
    cov_sample_proj = compute_raw_covariance(
        raw_sample.copy().apply_proj())

    cov_sss = compute_raw_covariance(raw_sss)
    cov_sss_proj = compute_raw_covariance(
        raw_sss.copy().apply_proj())

    picks_all_sample = pick_types(raw_sample.info, meg=True, eeg=True)
    picks_all_sss = pick_types(raw_sss.info, meg=True, eeg=True)

    info_sample = pick_info(raw_sample.info, picks_all_sample)
    picks_stack_sample = [('eeg', pick_types(info_sample, meg=False,
                                             eeg=True))]
    picks_stack_sample += [('meg', pick_types(info_sample, meg=True))]
    picks_stack_sample += [('all',
                            pick_types(info_sample, meg=True, eeg=True))]

    info_sss = pick_info(raw_sss.info, picks_all_sss)
    picks_stack_somato = [('eeg', pick_types(info_sss, meg=False, eeg=True))]
    picks_stack_somato += [('meg', pick_types(info_sss, meg=True))]
    picks_stack_somato += [('all',
                            pick_types(info_sss, meg=True, eeg=True))]

    iter_tests = list(itt.product(
        [(cov_sample, picks_stack_sample, info_sample),
         (cov_sample_proj, picks_stack_sample, info_sample),
         (cov_sss, picks_stack_somato, info_sss),
         (cov_sss_proj, picks_stack_somato, info_sss)],  # sss
        [dict(mag=1e15, grad=1e13, eeg=1e6)]
    ))

    for (cov, picks_list, this_info), scalings in iter_tests:
        for ch_type, picks in picks_list:

            this_very_info = pick_info(this_info, picks)

            # compute subset of projs
            this_projs = [c['active'] and
                          len(set(c['data']['col_names'])
                              .intersection(set(this_very_info['ch_names']))) >
                          0 for c in cov['projs']]
            n_projs = sum(this_projs)

            # count channel types
            ch_types = [channel_type(this_very_info, idx)
                        for idx in range(len(picks))]
            n_eeg, n_mag, n_grad = [ch_types.count(k) for k in
                                    ['eeg', 'mag', 'grad']]
            n_meg = n_mag + n_grad
            if ch_type in ('all', 'eeg'):
                n_projs_eeg = 1
            else:
                n_projs_eeg = 0

            # check sss
            if 'proc_history' in this_very_info:
                mf = this_very_info['proc_history'][0]['max_info']
                n_free = _get_sss_rank(mf)
                if 'mag' not in ch_types and 'grad' not in ch_types:
                    n_free = 0
                # - n_projs XXX clarify
                expected_rank = n_free + n_eeg
                if n_projs > 0 and ch_type in ('all', 'eeg'):
                    expected_rank -= n_projs_eeg
            else:
                expected_rank = n_meg + n_eeg - n_projs

            C = cov['data'][np.ix_(picks, picks)]
            est_rank = _estimate_rank_meeg_cov(C, this_very_info,
                                               scalings=scalings)

            assert_equal(expected_rank, est_rank)

Example 33

Project: AWS-Lambda-ML-Microservice-Skeleton Source File: test_index_tricks.py
    def test_1d_only(self):
        idx2d = [[1, 2, 3], [4, 5, 6]]
        assert_raises(ValueError, np.ix_, idx2d)