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
3
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
3
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
3
Example 3
@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
3
Example 4
@property
def V(self):
"""
Get the KxK matrix of scales
:return:
"""
return self.v[np.ix_(self.c, self.c)]
3
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)
3
Example 6
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)
3
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
3
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
3
Example 9
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))
0
Example 10
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
0
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
0
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
0
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,
),
}
0
Example 14
Project: qutip Source File: test_sparse.py
def _permutateIndexes(array, row_perm, col_perm):
return array[np.ix_(row_perm, col_perm)]
0
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
0
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
0
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
0
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
0
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
0
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()
0
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()
0
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()
0
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()
0
Example 24
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
0
Example 25
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
0
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)]
0
Example 27
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()
0
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
0
Example 29
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)]
0
Example 30
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
0
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
0
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)
0
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)