https://github.com/PyPSA/PyPSA
Tip revision: c2788daba46056552fd60f5bb1155342e22278a2 authored by Fabian on 13 October 2022, 14:10:19 UTC
doc: add a comment to optimal_power_flow section
doc: add a comment to optimal_power_flow section
Tip revision: c2788da
pf.py
# -*- coding: utf-8 -*-
"""
Power flow functionality.
"""
__author__ = (
"PyPSA Developers, see https://pypsa.readthedocs.io/en/latest/developers.html"
)
__copyright__ = (
"Copyright 2015-2022 PyPSA Developers, see https://pypsa.readthedocs.io/en/latest/developers.html, "
"MIT License"
)
import logging
logger = logging.getLogger(__name__)
import time
from operator import itemgetter
import networkx as nx
import numpy as np
import pandas as pd
from numpy import ones, r_
from numpy.linalg import norm
from pandas.api.types import is_list_like
from scipy.sparse import csc_matrix, csr_matrix, dok_matrix
from scipy.sparse import hstack as shstack
from scipy.sparse import issparse
from scipy.sparse import vstack as svstack
from scipy.sparse.linalg import spsolve
from pypsa.descriptors import (
Dict,
allocate_series_dataframes,
degree,
get_switchable_as_dense,
zsum,
)
pd.Series.zsum = zsum
def normed(s):
return s / s.sum()
def real(X):
return np.real(X.to_numpy())
def imag(X):
return np.imag(X.to_numpy())
def _as_snapshots(network, snapshots):
if snapshots is None:
snapshots = network.snapshots
if not is_list_like(snapshots):
snapshots = pd.Index([snapshots])
if not isinstance(snapshots, pd.MultiIndex):
snapshots = pd.Index(snapshots)
assert snapshots.isin(network.snapshots).all()
snapshots.name = "snapshot"
return snapshots
def _allocate_pf_outputs(network, linear=False):
to_allocate = {
"Generator": ["p"],
"Load": ["p"],
"StorageUnit": ["p"],
"Store": ["p"],
"ShuntImpedance": ["p"],
"Bus": ["p", "v_ang", "v_mag_pu"],
"Line": ["p0", "p1"],
"Transformer": ["p0", "p1"],
"Link": ["p" + col[3:] for col in network.links.columns if col[:3] == "bus"],
}
if not linear:
for component, attrs in to_allocate.items():
if "p" in attrs:
attrs.append("q")
if "p0" in attrs and component != "Link":
attrs.extend(["q0", "q1"])
allocate_series_dataframes(network, to_allocate)
def _calculate_controllable_nodal_power_balance(
sub_network, network, snapshots, buses_o
):
for n in ("q", "p"):
# allow all one ports to dispatch as set
for c in sub_network.iterate_components(
network.controllable_one_port_components
):
c_n_set = get_switchable_as_dense(
network, c.name, n + "_set", snapshots, c.ind
)
network.pnl(c.name)[n].loc[snapshots, c.ind] = c_n_set
# set the power injection at each node from controllable components
network.buses_t[n].loc[snapshots, buses_o] = sum(
[
(
(c.pnl[n].loc[snapshots, c.ind] * c.df.loc[c.ind, "sign"])
.groupby(c.df.loc[c.ind, "bus"], axis=1)
.sum()
.reindex(columns=buses_o, fill_value=0.0)
)
for c in sub_network.iterate_components(
network.controllable_one_port_components
)
]
)
if n == "p":
network.buses_t[n].loc[snapshots, buses_o] += sum(
[
(
-c.pnl[n + str(i)]
.loc[snapshots]
.groupby(c.df["bus" + str(i)], axis=1)
.sum()
.reindex(columns=buses_o, fill_value=0)
)
for c in network.iterate_components(
network.controllable_branch_components
)
for i in [int(col[3:]) for col in c.df.columns if col[:3] == "bus"]
]
)
def _network_prepare_and_run_pf(
network,
snapshots,
skip_pre,
linear=False,
distribute_slack=False,
slack_weights="p_set",
**kwargs
):
if linear:
sub_network_pf_fun = sub_network_lpf
sub_network_prepare_fun = calculate_B_H
else:
sub_network_pf_fun = sub_network_pf
sub_network_prepare_fun = calculate_Y
if not skip_pre:
network.determine_network_topology()
calculate_dependent_values(network)
_allocate_pf_outputs(network, linear)
snapshots = _as_snapshots(network, snapshots)
# deal with links
if not network.links.empty:
p_set = get_switchable_as_dense(network, "Link", "p_set", snapshots)
network.links_t.p0.loc[snapshots] = p_set.loc[snapshots]
for i in [
int(col[3:])
for col in network.links.columns
if col[:3] == "bus" and col != "bus0"
]:
eff_name = "efficiency" if i == 1 else "efficiency{}".format(i)
efficiency = get_switchable_as_dense(network, "Link", eff_name, snapshots)
links = network.links.index[network.links["bus{}".format(i)] != ""]
network.links_t["p{}".format(i)].loc[snapshots, links] = (
-network.links_t.p0.loc[snapshots, links]
* efficiency.loc[snapshots, links]
)
itdf = pd.DataFrame(index=snapshots, columns=network.sub_networks.index, dtype=int)
difdf = pd.DataFrame(index=snapshots, columns=network.sub_networks.index)
cnvdf = pd.DataFrame(
index=snapshots, columns=network.sub_networks.index, dtype=bool
)
for sub_network in network.sub_networks.obj:
if not skip_pre:
find_bus_controls(sub_network)
branches_i = sub_network.branches_i()
if len(branches_i) > 0:
sub_network_prepare_fun(sub_network, skip_pre=True)
if isinstance(slack_weights, dict):
sn_slack_weights = slack_weights[sub_network.name]
else:
sn_slack_weights = slack_weights
if isinstance(sn_slack_weights, dict):
sn_slack_weights = pd.Series(sn_slack_weights)
if not linear:
# escape for single-bus sub-network
if len(sub_network.buses()) <= 1:
(
itdf[sub_network.name],
difdf[sub_network.name],
cnvdf[sub_network.name],
) = sub_network_pf_singlebus(
sub_network,
snapshots=snapshots,
skip_pre=True,
distribute_slack=distribute_slack,
slack_weights=sn_slack_weights,
)
else:
(
itdf[sub_network.name],
difdf[sub_network.name],
cnvdf[sub_network.name],
) = sub_network_pf_fun(
sub_network,
snapshots=snapshots,
skip_pre=True,
distribute_slack=distribute_slack,
slack_weights=sn_slack_weights,
**kwargs
)
else:
sub_network_pf_fun(
sub_network, snapshots=snapshots, skip_pre=True, **kwargs
)
if not linear:
return Dict({"n_iter": itdf, "error": difdf, "converged": cnvdf})
def network_pf(
network,
snapshots=None,
skip_pre=False,
x_tol=1e-6,
use_seed=False,
distribute_slack=False,
slack_weights="p_set",
):
"""
Full non-linear power flow for generic network.
Parameters
----------
snapshots : list-like|single snapshot
A subset or an elements of network.snapshots on which to run
the power flow, defaults to network.snapshots
skip_pre : bool, default False
Skip the preliminary steps of computing topology, calculating dependent values and finding bus controls.
x_tol: float
Tolerance for Newton-Raphson power flow.
use_seed : bool, default False
Use a seed for the initial guess for the Newton-Raphson algorithm.
distribute_slack : bool, default False
If ``True``, distribute the slack power across generators proportional to generator dispatch by default
or according to the distribution scheme provided in ``slack_weights``.
If ``False`` only the slack generator takes up the slack.
slack_weights : dict|str, default 'p_set'
Distribution scheme describing how to determine the fraction of the total slack power
(of each sub network individually) a bus of the subnetwork takes up.
Default is to distribute proportional to generator dispatch ('p_set').
Another option is to distribute proportional to (optimised) nominal capacity ('p_nom' or 'p_nom_opt').
Custom weights can be specified via a dictionary that has a key for each
subnetwork index (``network.sub_networks.index``) and a
pandas.Series/dict with buses or generators of the
corresponding subnetwork as index/keys.
When specifying custom weights with buses as index/keys the slack power of a bus is distributed
among its generators in proportion to their nominal capacity (``p_nom``) if given, otherwise evenly.
Returns
-------
dict
Dictionary with keys 'n_iter', 'converged', 'error' and dataframe
values indicating number of iterations, convergence status, and
iteration error for each snapshot (rows) and sub_network (columns)
"""
return _network_prepare_and_run_pf(
network,
snapshots,
skip_pre,
linear=False,
x_tol=x_tol,
use_seed=use_seed,
distribute_slack=distribute_slack,
slack_weights=slack_weights,
)
def newton_raphson_sparse(
f,
guess,
dfdx,
x_tol=1e-10,
lim_iter=100,
distribute_slack=False,
slack_weights=None,
):
"""Solve f(x) = 0 with initial guess for x and dfdx(x). dfdx(x) should
return a sparse Jacobian. Terminate if error on norm of f(x) is <
x_tol or there were more than lim_iter iterations.
"""
slack_args = {"distribute_slack": distribute_slack, "slack_weights": slack_weights}
converged = False
n_iter = 0
F = f(guess, **slack_args)
diff = norm(F, np.Inf)
logger.debug("Error at iteration %d: %f", n_iter, diff)
while diff > x_tol and n_iter < lim_iter:
n_iter += 1
guess = guess - spsolve(dfdx(guess, **slack_args), F)
F = f(guess, **slack_args)
diff = norm(F, np.Inf)
logger.debug("Error at iteration %d: %f", n_iter, diff)
if diff > x_tol:
logger.warning(
'Warning, we didn\'t reach the required tolerance within %d iterations, error is at %f. See the section "Troubleshooting" in the documentation for tips to fix this. ',
n_iter,
diff,
)
elif not np.isnan(diff):
converged = True
return guess, n_iter, diff, converged
def sub_network_pf_singlebus(
sub_network,
snapshots=None,
skip_pre=False,
distribute_slack=False,
slack_weights="p_set",
linear=False,
):
"""
Non-linear power flow for a sub-network consiting of a single bus.
Parameters
----------
snapshots : list-like|single snapshot
A subset or an elements of network.snapshots on which to run
the power flow, defaults to network.snapshots
skip_pre: bool, default False
Skip the preliminary steps of computing topology, calculating dependent values and finding bus controls.
distribute_slack : bool, default False
If ``True``, distribute the slack power across generators proportional to generator dispatch by default
or according to the distribution scheme provided in ``slack_weights``.
If ``False`` only the slack generator takes up the slack.
slack_weights : pandas.Series|str, default 'p_set'
Distribution scheme describing how to determine the fraction of the total slack power
a bus of the subnetwork takes up. Default is to distribute proportional to generator dispatch
('p_set'). Another option is to distribute proportional to (optimised) nominal capacity ('p_nom' or 'p_nom_opt').
Custom weights can be provided via a pandas.Series/dict
that has the generators of the single bus as index/keys.
"""
snapshots = _as_snapshots(sub_network.network, snapshots)
network = sub_network.network
logger.info(
"Balancing power on single-bus sub-network {} for snapshots {}".format(
sub_network, snapshots
)
)
if not skip_pre:
find_bus_controls(sub_network)
_allocate_pf_outputs(network, linear=False)
if isinstance(slack_weights, dict):
slack_weights = pd.Series(slack_weights)
buses_o = sub_network.buses_o
_calculate_controllable_nodal_power_balance(
sub_network, network, snapshots, buses_o
)
v_mag_pu_set = get_switchable_as_dense(network, "Bus", "v_mag_pu_set", snapshots)
network.buses_t.v_mag_pu.loc[snapshots, sub_network.slack_bus] = v_mag_pu_set.loc[
:, sub_network.slack_bus
]
network.buses_t.v_ang.loc[snapshots, sub_network.slack_bus] = 0.0
if distribute_slack:
for bus, group in sub_network.generators().groupby("bus"):
if slack_weights in ["p_nom", "p_nom_opt"]:
assert (
not all(network.generators[slack_weights]) == 0
), "Invalid slack weights! Generator attribute {} is always zero.".format(
slack_weights
)
bus_generator_shares = (
network.generators[slack_weights]
.loc[group.index]
.pipe(normed)
.fillna(0)
)
elif slack_weights == "p_set":
generators_t_p_choice = get_switchable_as_dense(
network, "Generator", slack_weights, snapshots
)
assert (
not generators_t_p_choice.isna().all().all()
), "Invalid slack weights! Generator attribute {} is always NaN.".format(
slack_weights
)
assert (
not (generators_t_p_choice == 0).all().all()
), "Invalid slack weights! Generator attribute {} is always zero.".format(
slack_weights
)
bus_generator_shares = (
generators_t_p_choice.loc[snapshots, group.index]
.apply(normed, axis=1)
.fillna(0)
)
else:
bus_generator_shares = slack_weights.pipe(normed).fillna(0)
network.generators_t.p.loc[
snapshots, group.index
] += bus_generator_shares.multiply(
-network.buses_t.p.loc[snapshots, bus], axis=0
)
else:
network.generators_t.p.loc[
snapshots, sub_network.slack_generator
] -= network.buses_t.p.loc[snapshots, sub_network.slack_bus]
network.generators_t.q.loc[
snapshots, sub_network.slack_generator
] -= network.buses_t.q.loc[snapshots, sub_network.slack_bus]
network.buses_t.p.loc[snapshots, sub_network.slack_bus] = 0.0
network.buses_t.q.loc[snapshots, sub_network.slack_bus] = 0.0
return 0, 0.0, True # dummy substitute for newton raphson output
def sub_network_pf(
sub_network,
snapshots=None,
skip_pre=False,
x_tol=1e-6,
use_seed=False,
distribute_slack=False,
slack_weights="p_set",
):
"""
Non-linear power flow for connected sub-network.
Parameters
----------
snapshots : list-like|single snapshot
A subset or an elements of network.snapshots on which to run
the power flow, defaults to network.snapshots
skip_pre: bool, default False
Skip the preliminary steps of computing topology, calculating dependent values and finding bus controls.
x_tol: float
Tolerance for Newton-Raphson power flow.
use_seed : bool, default False
Use a seed for the initial guess for the Newton-Raphson algorithm.
distribute_slack : bool, default False
If ``True``, distribute the slack power across generators proportional to generator dispatch by default
or according to the distribution scheme provided in ``slack_weights``.
If ``False`` only the slack generator takes up the slack.
slack_weights : pandas.Series|str, default 'p_set'
Distribution scheme describing how to determine the fraction of the total slack power
a bus of the subnetwork takes up. Default is to distribute proportional to generator dispatch
('p_set'). Another option is to distribute proportional to (optimised) nominal capacity ('p_nom' or 'p_nom_opt').
Custom weights can be provided via a pandas.Series/dict
that has the buses or the generators of the subnetwork as index/keys.
When using custom weights with buses as index/keys the slack power of a bus is distributed
among its generators in proportion to their nominal capacity (``p_nom``) if given, otherwise evenly.
Returns
-------
Tuple of three pandas.Series indicating number of iterations,
remaining error, and convergence status for each snapshot
"""
assert isinstance(
slack_weights, (str, pd.Series, dict)
), "Type of 'slack_weights' must be string, pd.Series or dict. Is {}.".format(
type(slack_weights)
)
if isinstance(slack_weights, dict):
slack_weights = pd.Series(slack_weights)
elif isinstance(slack_weights, str):
valid_strings = ["p_nom", "p_nom_opt", "p_set"]
assert (
slack_weights in valid_strings
), "String value for 'slack_weights' must be one of {}. Is {}.".format(
valid_strings, slack_weights
)
snapshots = _as_snapshots(sub_network.network, snapshots)
logger.info(
"Performing non-linear load-flow on {} sub-network {} for snapshots {}".format(
sub_network.network.sub_networks.at[sub_network.name, "carrier"],
sub_network,
snapshots,
)
)
# _sub_network_prepare_pf(sub_network, snapshots, skip_pre, calculate_Y)
network = sub_network.network
if not skip_pre:
calculate_dependent_values(network)
find_bus_controls(sub_network)
_allocate_pf_outputs(network, linear=False)
# get indices for the components on this subnetwork
branches_i = sub_network.branches_i()
buses_o = sub_network.buses_o
sn_buses = sub_network.buses().index
sn_generators = sub_network.generators().index
generator_slack_weights_b = False
bus_slack_weights_b = False
if isinstance(slack_weights, pd.Series):
if all(i in sn_generators for i in slack_weights.index):
generator_slack_weights_b = True
elif all(i in sn_buses for i in slack_weights.index):
bus_slack_weights_b = True
else:
raise AssertionError(
"Custom slack weights pd.Series/dict must only have the",
"generators or buses of the subnetwork as index/keys.",
)
if not skip_pre and len(branches_i) > 0:
calculate_Y(sub_network, skip_pre=True)
_calculate_controllable_nodal_power_balance(
sub_network, network, snapshots, buses_o
)
def f(guess, distribute_slack=False, slack_weights=None):
last_pq = -1 if distribute_slack else None
network.buses_t.v_ang.loc[now, sub_network.pvpqs] = guess[
: len(sub_network.pvpqs)
]
network.buses_t.v_mag_pu.loc[now, sub_network.pqs] = guess[
len(sub_network.pvpqs) : last_pq
]
v_mag_pu = network.buses_t.v_mag_pu.loc[now, buses_o]
v_ang = network.buses_t.v_ang.loc[now, buses_o]
V = v_mag_pu * np.exp(1j * v_ang)
if distribute_slack:
slack_power = slack_weights * guess[-1]
mismatch = V * np.conj(sub_network.Y * V) - s + slack_power
else:
mismatch = V * np.conj(sub_network.Y * V) - s
if distribute_slack:
F = r_[real(mismatch)[:], imag(mismatch)[1 + len(sub_network.pvs) :]]
else:
F = r_[real(mismatch)[1:], imag(mismatch)[1 + len(sub_network.pvs) :]]
return F
def dfdx(guess, distribute_slack=False, slack_weights=None):
last_pq = -1 if distribute_slack else None
network.buses_t.v_ang.loc[now, sub_network.pvpqs] = guess[
: len(sub_network.pvpqs)
]
network.buses_t.v_mag_pu.loc[now, sub_network.pqs] = guess[
len(sub_network.pvpqs) : last_pq
]
v_mag_pu = network.buses_t.v_mag_pu.loc[now, buses_o]
v_ang = network.buses_t.v_ang.loc[now, buses_o]
V = v_mag_pu * np.exp(1j * v_ang)
index = r_[: len(buses_o)]
# make sparse diagonal matrices
V_diag = csr_matrix((V, (index, index)))
V_norm_diag = csr_matrix((V / abs(V), (index, index)))
I_diag = csr_matrix((sub_network.Y * V, (index, index)))
dS_dVa = 1j * V_diag * np.conj(I_diag - sub_network.Y * V_diag)
dS_dVm = V_norm_diag * np.conj(I_diag) + V_diag * np.conj(
sub_network.Y * V_norm_diag
)
J10 = dS_dVa[1 + len(sub_network.pvs) :, 1:].imag
J11 = dS_dVm[1 + len(sub_network.pvs) :, 1 + len(sub_network.pvs) :].imag
if distribute_slack:
J00 = dS_dVa[:, 1:].real
J01 = dS_dVm[:, 1 + len(sub_network.pvs) :].real
J02 = csr_matrix(slack_weights, (1, 1 + len(sub_network.pvpqs))).T
J12 = csr_matrix((1, len(sub_network.pqs))).T
J_P_blocks = [J00, J01, J02]
J_Q_blocks = [J10, J11, J12]
else:
J00 = dS_dVa[1:, 1:].real
J01 = dS_dVm[1:, 1 + len(sub_network.pvs) :].real
J_P_blocks = [J00, J01]
J_Q_blocks = [J10, J11]
J = svstack([shstack(J_P_blocks), shstack(J_Q_blocks)], format="csr")
return J
# Set what we know: slack V and v_mag_pu for PV buses
v_mag_pu_set = get_switchable_as_dense(network, "Bus", "v_mag_pu_set", snapshots)
network.buses_t.v_mag_pu.loc[snapshots, sub_network.pvs] = v_mag_pu_set.loc[
:, sub_network.pvs
]
network.buses_t.v_mag_pu.loc[snapshots, sub_network.slack_bus] = v_mag_pu_set.loc[
:, sub_network.slack_bus
]
network.buses_t.v_ang.loc[snapshots, sub_network.slack_bus] = 0.0
if not use_seed:
network.buses_t.v_mag_pu.loc[snapshots, sub_network.pqs] = 1.0
network.buses_t.v_ang.loc[snapshots, sub_network.pvpqs] = 0.0
slack_args = {"distribute_slack": distribute_slack}
slack_variable_b = 1 if distribute_slack else 0
if distribute_slack:
if isinstance(slack_weights, str) and slack_weights == "p_set":
generators_t_p_choice = get_switchable_as_dense(
network, "Generator", slack_weights, snapshots
)
bus_generation = generators_t_p_choice.rename(
columns=network.generators.bus
)
slack_weights_calc = (
pd.DataFrame(
bus_generation.groupby(bus_generation.columns, axis=1).sum(),
columns=buses_o,
)
.apply(normed, axis=1)
.fillna(0)
)
elif isinstance(slack_weights, str) and slack_weights in ["p_nom", "p_nom_opt"]:
assert (
not all(network.generators[slack_weights]) == 0
), "Invalid slack weights! Generator attribute {} is always zero.".format(
slack_weights
)
slack_weights_calc = (
network.generators.groupby("bus")
.sum()[slack_weights]
.reindex(buses_o)
.pipe(normed)
.fillna(0)
)
elif generator_slack_weights_b:
# convert generator-based slack weights to bus-based slack weights
slack_weights_calc = (
slack_weights.rename(network.generators.bus)
.groupby(slack_weights.index.name)
.sum()
.reindex(buses_o)
.pipe(normed)
.fillna(0)
)
elif bus_slack_weights_b:
# take bus-based slack weights
slack_weights_calc = slack_weights.reindex(buses_o).pipe(normed).fillna(0)
ss = np.empty((len(snapshots), len(buses_o)), dtype=complex)
roots = np.empty(
(
len(snapshots),
len(sub_network.pvpqs) + len(sub_network.pqs) + slack_variable_b,
)
)
iters = pd.Series(0, index=snapshots)
diffs = pd.Series(index=snapshots, dtype=float)
convs = pd.Series(False, index=snapshots)
for i, now in enumerate(snapshots):
p = network.buses_t.p.loc[now, buses_o]
q = network.buses_t.q.loc[now, buses_o]
ss[i] = s = p + 1j * q
# Make a guess for what we don't know: V_ang for PV and PQs and v_mag_pu for PQ buses
guess = r_[
network.buses_t.v_ang.loc[now, sub_network.pvpqs],
network.buses_t.v_mag_pu.loc[now, sub_network.pqs],
]
if distribute_slack:
guess = np.append(guess, [0]) # for total slack power
if isinstance(slack_weights, str) and slack_weights == "p_set":
# snapshot-dependent slack weights
slack_args["slack_weights"] = slack_weights_calc.loc[now]
else:
slack_args["slack_weights"] = slack_weights_calc
# Now try and solve
start = time.time()
roots[i], n_iter, diff, converged = newton_raphson_sparse(
f, guess, dfdx, x_tol=x_tol, **slack_args
)
logger.info(
"Newton-Raphson solved in %d iterations with error of %f in %f seconds",
n_iter,
diff,
time.time() - start,
)
iters[now] = n_iter
diffs[now] = diff
convs[now] = converged
# now set everything
if distribute_slack:
last_pq = -1
slack_power = roots[:, -1]
else:
last_pq = None
network.buses_t.v_ang.loc[snapshots, sub_network.pvpqs] = roots[
:, : len(sub_network.pvpqs)
]
network.buses_t.v_mag_pu.loc[snapshots, sub_network.pqs] = roots[
:, len(sub_network.pvpqs) : last_pq
]
v_mag_pu = network.buses_t.v_mag_pu.loc[snapshots, buses_o].values
v_ang = network.buses_t.v_ang.loc[snapshots, buses_o].values
V = v_mag_pu * np.exp(1j * v_ang)
# add voltages to branches
buses_indexer = buses_o.get_indexer
branch_bus0 = []
branch_bus1 = []
for c in sub_network.iterate_components(network.passive_branch_components):
branch_bus0 += list(c.df.loc[c.ind, "bus0"])
branch_bus1 += list(c.df.loc[c.ind, "bus1"])
v0 = V[:, buses_indexer(branch_bus0)]
v1 = V[:, buses_indexer(branch_bus1)]
i0 = np.empty((len(snapshots), sub_network.Y0.shape[0]), dtype=complex)
i1 = np.empty((len(snapshots), sub_network.Y1.shape[0]), dtype=complex)
for i, now in enumerate(snapshots):
i0[i] = sub_network.Y0 * V[i]
i1[i] = sub_network.Y1 * V[i]
s0 = pd.DataFrame(v0 * np.conj(i0), columns=branches_i, index=snapshots)
s1 = pd.DataFrame(v1 * np.conj(i1), columns=branches_i, index=snapshots)
for c in sub_network.iterate_components(network.passive_branch_components):
s0t = s0.loc[:, c.name]
s1t = s1.loc[:, c.name]
network.pnl(c.name).p0.loc[snapshots, s0t.columns] = s0t.values.real
network.pnl(c.name).q0.loc[snapshots, s0t.columns] = s0t.values.imag
network.pnl(c.name).p1.loc[snapshots, s1t.columns] = s1t.values.real
network.pnl(c.name).q1.loc[snapshots, s1t.columns] = s1t.values.imag
s_calc = np.empty((len(snapshots), len(buses_o)), dtype=complex)
for i in np.arange(len(snapshots)):
s_calc[i] = V[i] * np.conj(sub_network.Y * V[i])
slack_index = buses_o.get_loc(sub_network.slack_bus)
if distribute_slack:
network.buses_t.p.loc[snapshots, sn_buses] = s_calc.real[
:, buses_indexer(sn_buses)
]
else:
network.buses_t.p.loc[snapshots, sub_network.slack_bus] = s_calc[
:, slack_index
].real
network.buses_t.q.loc[snapshots, sub_network.slack_bus] = s_calc[
:, slack_index
].imag
network.buses_t.q.loc[snapshots, sub_network.pvs] = s_calc[
:, buses_indexer(sub_network.pvs)
].imag
# set shunt impedance powers
shunt_impedances_i = sub_network.shunt_impedances_i()
if len(shunt_impedances_i):
# add voltages
shunt_impedances_v_mag_pu = v_mag_pu[
:, buses_indexer(network.shunt_impedances.loc[shunt_impedances_i, "bus"])
]
network.shunt_impedances_t.p.loc[snapshots, shunt_impedances_i] = (
shunt_impedances_v_mag_pu**2
) * network.shunt_impedances.loc[shunt_impedances_i, "g_pu"].values
network.shunt_impedances_t.q.loc[snapshots, shunt_impedances_i] = (
shunt_impedances_v_mag_pu**2
) * network.shunt_impedances.loc[shunt_impedances_i, "b_pu"].values
# let slack generator take up the slack
if distribute_slack:
distributed_slack_power = (
network.buses_t.p.loc[snapshots, sn_buses]
- ss[:, buses_indexer(sn_buses)].real
)
for bus, group in sub_network.generators().groupby("bus"):
if isinstance(slack_weights, str) and slack_weights == "p_set":
generators_t_p_choice = get_switchable_as_dense(
network, "Generator", slack_weights, snapshots
)
bus_generator_shares = (
generators_t_p_choice.loc[snapshots, group.index]
.apply(normed, axis=1)
.fillna(0)
)
network.generators_t.p.loc[
snapshots, group.index
] += bus_generator_shares.multiply(
distributed_slack_power.loc[snapshots, bus], axis=0
)
else:
if generator_slack_weights_b:
bus_generator_shares = (
slack_weights.loc[group.index].pipe(normed).fillna(0)
)
else:
bus_generators_p_nom = network.generators.p_nom.loc[group.index]
# distribute evenly if no p_nom given
if all(bus_generators_p_nom) == 0:
bus_generators_p_nom = 1
bus_generator_shares = bus_generators_p_nom.pipe(normed).fillna(0)
network.generators_t.p.loc[
snapshots, group.index
] += distributed_slack_power.loc[snapshots, bus].apply(
lambda row: row * bus_generator_shares
)
else:
network.generators_t.p.loc[snapshots, sub_network.slack_generator] += (
network.buses_t.p.loc[snapshots, sub_network.slack_bus]
- ss[:, slack_index].real
)
# set the Q of the slack and PV generators
network.generators_t.q.loc[snapshots, sub_network.slack_generator] += (
network.buses_t.q.loc[snapshots, sub_network.slack_bus]
- ss[:, slack_index].imag
)
network.generators_t.q.loc[
snapshots, network.buses.loc[sub_network.pvs, "generator"]
] += np.asarray(
network.buses_t.q.loc[snapshots, sub_network.pvs]
- ss[:, buses_indexer(sub_network.pvs)].imag
)
return iters, diffs, convs
def network_lpf(network, snapshots=None, skip_pre=False):
"""
Linear power flow for generic network.
Parameters
----------
snapshots : list-like|single snapshot
A subset or an elements of network.snapshots on which to run
the power flow, defaults to network.snapshots
skip_pre : bool, default False
Skip the preliminary steps of computing topology, calculating
dependent values and finding bus controls.
Returns
-------
None
"""
_network_prepare_and_run_pf(network, snapshots, skip_pre, linear=True)
def apply_line_types(network):
"""
Calculate line electrical parameters x, r, b, g from standard types.
"""
lines_with_types_b = network.lines.type != ""
if lines_with_types_b.zsum() == 0:
return
missing_types = pd.Index(
network.lines.loc[lines_with_types_b, "type"].unique()
).difference(network.line_types.index)
assert (
missing_types.empty
), "The type(s) {} do(es) not exist in network.line_types".format(
", ".join(missing_types)
)
# Get a copy of the lines data
l = network.lines.loc[lines_with_types_b, ["type", "length", "num_parallel"]].join(
network.line_types, on="type"
)
for attr in ["r", "x"]:
l[attr] = l[attr + "_per_length"] * l["length"] / l["num_parallel"]
l["b"] = (
2
* np.pi
* 1e-9
* l["f_nom"]
* l["c_per_length"]
* l["length"]
* l["num_parallel"]
)
# now set calculated values on live lines
for attr in ["r", "x", "b"]:
network.lines.loc[lines_with_types_b, attr] = l[attr]
def apply_transformer_types(network):
"""
Calculate transformer electrical parameters x, r, b, g from standard types.
"""
trafos_with_types_b = network.transformers.type != ""
if trafos_with_types_b.zsum() == 0:
return
missing_types = pd.Index(
network.transformers.loc[trafos_with_types_b, "type"].unique()
).difference(network.transformer_types.index)
assert (
missing_types.empty
), "The type(s) {} do(es) not exist in network.transformer_types".format(
", ".join(missing_types)
)
# Get a copy of the transformers data
# (joining pulls in "phase_shift", "s_nom", "tap_side" from TransformerType)
t = network.transformers.loc[
trafos_with_types_b, ["type", "tap_position", "num_parallel"]
].join(network.transformer_types, on="type")
t["r"] = t["vscr"] / 100.0
t["x"] = np.sqrt((t["vsc"] / 100.0) ** 2 - t["r"] ** 2)
# NB: b and g are per unit of s_nom
t["g"] = t["pfe"] / (1000.0 * t["s_nom"])
# for some bizarre reason, some of the standard types in pandapower have i0^2 < g^2
t["b"] = -np.sqrt(((t["i0"] / 100.0) ** 2 - t["g"] ** 2).clip(lower=0))
for attr in ["r", "x"]:
t[attr] /= t["num_parallel"]
for attr in ["b", "g"]:
t[attr] *= t["num_parallel"]
# deal with tap positions
t["tap_ratio"] = 1.0 + (t["tap_position"] - t["tap_neutral"]) * (
t["tap_step"] / 100.0
)
# now set calculated values on live transformers
for attr in ["r", "x", "g", "b", "phase_shift", "s_nom", "tap_side", "tap_ratio"]:
network.transformers.loc[trafos_with_types_b, attr] = t[attr]
# TODO: status, rate_A
def wye_to_delta(z1, z2, z3):
"""
Follows http://home.earthlink.net/~w6rmk/math/wyedelta.htm.
"""
summand = z1 * z2 + z2 * z3 + z3 * z1
return (summand / z2, summand / z1, summand / z3)
def apply_transformer_t_model(network):
"""
Convert given T-model parameters to PI-model parameters using wye-delta
transformation.
"""
z_series = network.transformers.r_pu + 1j * network.transformers.x_pu
y_shunt = network.transformers.g_pu + 1j * network.transformers.b_pu
ts_b = (network.transformers.model == "t") & (y_shunt != 0.0)
if ts_b.zsum() == 0:
return
za, zb, zc = wye_to_delta(
z_series.loc[ts_b] / 2, z_series.loc[ts_b] / 2, 1 / y_shunt.loc[ts_b]
)
network.transformers.loc[ts_b, "r_pu"] = real(zc)
network.transformers.loc[ts_b, "x_pu"] = imag(zc)
network.transformers.loc[ts_b, "g_pu"] = real(2 / za)
network.transformers.loc[ts_b, "b_pu"] = imag(2 / za)
def calculate_dependent_values(network):
"""
Calculate per unit impedances and append voltages to lines and shunt
impedances.
"""
apply_line_types(network)
apply_transformer_types(network)
network.lines["v_nom"] = network.lines.bus0.map(network.buses.v_nom)
network.lines.loc[network.lines.carrier == "", "carrier"] = network.lines.bus0.map(
network.buses.carrier
)
network.lines["x_pu"] = network.lines.x / (network.lines.v_nom**2)
network.lines["r_pu"] = network.lines.r / (network.lines.v_nom**2)
network.lines["b_pu"] = network.lines.b * network.lines.v_nom**2
network.lines["g_pu"] = network.lines.g * network.lines.v_nom**2
network.lines["x_pu_eff"] = network.lines["x_pu"]
network.lines["r_pu_eff"] = network.lines["r_pu"]
# convert transformer impedances from base power s_nom to base = 1 MVA
network.transformers["x_pu"] = network.transformers.x / network.transformers.s_nom
network.transformers["r_pu"] = network.transformers.r / network.transformers.s_nom
network.transformers["b_pu"] = network.transformers.b * network.transformers.s_nom
network.transformers["g_pu"] = network.transformers.g * network.transformers.s_nom
network.transformers["x_pu_eff"] = (
network.transformers["x_pu"] * network.transformers["tap_ratio"]
)
network.transformers["r_pu_eff"] = (
network.transformers["r_pu"] * network.transformers["tap_ratio"]
)
apply_transformer_t_model(network)
network.shunt_impedances["v_nom"] = network.shunt_impedances["bus"].map(
network.buses.v_nom
)
network.shunt_impedances["b_pu"] = (
network.shunt_impedances.b * network.shunt_impedances.v_nom**2
)
network.shunt_impedances["g_pu"] = (
network.shunt_impedances.g * network.shunt_impedances.v_nom**2
)
network.links.loc[network.links.carrier == "", "carrier"] = network.links.bus0.map(
network.buses.carrier
)
network.stores.loc[
network.stores.carrier == "", "carrier"
] = network.stores.bus.map(network.buses.carrier)
def find_slack_bus(sub_network):
"""
Find the slack bus in a connected sub-network.
"""
gens = sub_network.generators()
if len(gens) == 0:
# logger.warning("No generators in sub-network {}, better hope power is already balanced".format(sub_network.name))
sub_network.slack_generator = None
sub_network.slack_bus = sub_network.buses_i()[0]
else:
slacks = gens[gens.control == "Slack"].index
if len(slacks) == 0:
sub_network.slack_generator = gens.index[0]
sub_network.network.generators.loc[
sub_network.slack_generator, "control"
] = "Slack"
logger.debug(
"No slack generator found in sub-network {}, using {} as the slack generator".format(
sub_network.name, sub_network.slack_generator
)
)
elif len(slacks) == 1:
sub_network.slack_generator = slacks[0]
else:
sub_network.slack_generator = slacks[0]
sub_network.network.generators.loc[slacks[1:], "control"] = "PV"
logger.debug(
"More than one slack generator found in sub-network {}, using {} as the slack generator".format(
sub_network.name, sub_network.slack_generator
)
)
sub_network.slack_bus = gens.bus[sub_network.slack_generator]
# also put it into the dataframe
sub_network.network.sub_networks.at[
sub_network.name, "slack_bus"
] = sub_network.slack_bus
logger.debug(
"Slack bus for sub-network {} is {}".format(
sub_network.name, sub_network.slack_bus
)
)
def find_bus_controls(sub_network):
"""
Find slack and all PV and PQ buses for a sub_network.
This function also fixes sub_network.buses_o, a DataFrame ordered by
control type.
"""
network = sub_network.network
find_slack_bus(sub_network)
gens = sub_network.generators()
buses_i = sub_network.buses_i()
# default bus control is PQ
network.buses.loc[buses_i, "control"] = "PQ"
# find all buses with one or more gens with PV
pvs = gens[gens.control == "PV"].index.to_series()
if len(pvs) > 0:
pvs = pvs.groupby(gens.bus).first()
network.buses.loc[pvs.index, "control"] = "PV"
network.buses.loc[pvs.index, "generator"] = pvs
network.buses.loc[sub_network.slack_bus, "control"] = "Slack"
network.buses.loc[sub_network.slack_bus, "generator"] = sub_network.slack_generator
buses_control = network.buses.loc[buses_i, "control"]
sub_network.pvs = buses_control.index[buses_control == "PV"]
sub_network.pqs = buses_control.index[buses_control == "PQ"]
sub_network.pvpqs = sub_network.pvs.append(sub_network.pqs)
# order buses
sub_network.buses_o = sub_network.pvpqs.insert(0, sub_network.slack_bus)
def calculate_B_H(sub_network, skip_pre=False):
"""
Calculate B and H matrices for AC or DC sub-networks.
"""
network = sub_network.network
if not skip_pre:
calculate_dependent_values(network)
find_bus_controls(sub_network)
if network.sub_networks.at[sub_network.name, "carrier"] == "DC":
attribute = "r_pu_eff"
else:
attribute = "x_pu_eff"
# following leans heavily on pypower.makeBdc
# susceptances
b = 1.0 / np.concatenate(
[
(c.df.loc[c.ind, attribute]).values
for c in sub_network.iterate_components(network.passive_branch_components)
]
)
if np.isnan(b).any():
logger.warning(
"Warning! Some series impedances are zero - this will cause a singularity in LPF!"
)
b_diag = csr_matrix((b, (r_[: len(b)], r_[: len(b)])))
# incidence matrix
sub_network.K = sub_network.incidence_matrix(busorder=sub_network.buses_o)
sub_network.H = b_diag * sub_network.K.T
# weighted Laplacian
sub_network.B = sub_network.K * sub_network.H
sub_network.p_branch_shift = -b * np.concatenate(
[
(c.df.loc[c.ind, "phase_shift"]).values * np.pi / 180.0
if c.name == "Transformer"
else np.zeros((len(c.ind),))
for c in sub_network.iterate_components(network.passive_branch_components)
]
)
sub_network.p_bus_shift = sub_network.K * sub_network.p_branch_shift
def calculate_PTDF(sub_network, skip_pre=False):
"""
Calculate the Power Transfer Distribution Factor (PTDF) for sub_network.
Sets sub_network.PTDF as a (dense) numpy array.
Parameters
----------
sub_network : pypsa.SubNetwork
skip_pre : bool, default False
Skip the preliminary steps of computing topology, calculating dependent values,
finding bus controls and computing B and H.
"""
if not skip_pre:
calculate_B_H(sub_network)
# calculate inverse of B with slack removed
n_pvpq = len(sub_network.pvpqs)
index = np.r_[:n_pvpq]
I = csc_matrix((np.ones(n_pvpq), (index, index)))
B_inverse = spsolve(csc_matrix(sub_network.B[1:, 1:]), I)
# exception for two-node networks, where B_inverse is a 1d array
if issparse(B_inverse):
B_inverse = B_inverse.toarray()
elif B_inverse.shape == (1,):
B_inverse = B_inverse.reshape((1, 1))
# add back in zeroes for slack
B_inverse = np.hstack((np.zeros((n_pvpq, 1)), B_inverse))
B_inverse = np.vstack((np.zeros(n_pvpq + 1), B_inverse))
sub_network.PTDF = sub_network.H * B_inverse
def calculate_Y(sub_network, skip_pre=False):
"""
Calculate bus admittance matrices for AC sub-networks.
"""
if not skip_pre:
calculate_dependent_values(sub_network.network)
if sub_network.network.sub_networks.at[sub_network.name, "carrier"] != "AC":
logger.warning("Non-AC networks not supported for Y!")
return
branches = sub_network.branches()
buses_o = sub_network.buses_o
network = sub_network.network
# following leans heavily on pypower.makeYbus
# Copyright Richard Lincoln, Ray Zimmerman, BSD-style licence
num_branches = len(branches)
num_buses = len(buses_o)
y_se = 1 / (branches["r_pu"] + 1.0j * branches["x_pu"])
y_sh = branches["g_pu"] + 1.0j * branches["b_pu"]
tau = branches["tap_ratio"].fillna(1.0)
# catch some transformers falsely set with tau = 0 by pypower
tau[tau == 0] = 1.0
# define the HV tap ratios
tau_hv = pd.Series(1.0, branches.index)
tau_hv[branches.tap_side == 0] = tau[branches.tap_side == 0]
# define the LV tap ratios
tau_lv = pd.Series(1.0, branches.index)
tau_lv[branches.tap_side == 1] = tau[branches.tap_side == 1]
phase_shift = np.exp(1.0j * branches["phase_shift"].fillna(0.0) * np.pi / 180.0)
# build the admittance matrix elements for each branch
Y11 = (y_se + 0.5 * y_sh) / tau_lv**2
Y10 = -y_se / tau_lv / tau_hv / phase_shift
Y01 = -y_se / tau_lv / tau_hv / np.conj(phase_shift)
Y00 = (y_se + 0.5 * y_sh) / tau_hv**2
# bus shunt impedances
b_sh = (
network.shunt_impedances.b_pu.groupby(network.shunt_impedances.bus)
.sum()
.reindex(buses_o, fill_value=0.0)
)
g_sh = (
network.shunt_impedances.g_pu.groupby(network.shunt_impedances.bus)
.sum()
.reindex(buses_o, fill_value=0.0)
)
Y_sh = g_sh + 1.0j * b_sh
# get bus indices
bus0 = buses_o.get_indexer(branches.bus0)
bus1 = buses_o.get_indexer(branches.bus1)
# connection matrices
C0 = csr_matrix(
(ones(num_branches), (np.arange(num_branches), bus0)), (num_branches, num_buses)
)
C1 = csr_matrix(
(ones(num_branches), (np.arange(num_branches), bus1)), (num_branches, num_buses)
)
# build Y{0,1} such that Y{0,1} * V is the vector complex branch currents
i = r_[np.arange(num_branches), np.arange(num_branches)]
sub_network.Y0 = csr_matrix(
(r_[Y00, Y01], (i, r_[bus0, bus1])), (num_branches, num_buses)
)
sub_network.Y1 = csr_matrix(
(r_[Y10, Y11], (i, r_[bus0, bus1])), (num_branches, num_buses)
)
# now build bus admittance matrix
sub_network.Y = (
C0.T * sub_network.Y0
+ C1.T * sub_network.Y1
+ csr_matrix((Y_sh, (np.arange(num_buses), np.arange(num_buses))))
)
def aggregate_multi_graph(sub_network):
"""
Aggregate branches between same buses and replace with a single branch with
aggregated properties (e.g. s_nom is summed, length is averaged).
"""
network = sub_network.network
count = 0
seen = []
graph = sub_network.graph()
for u, v in graph.edges():
if (u, v) in seen:
continue
line_objs = list(graph.adj[u][v].keys())
if len(line_objs) > 1:
lines = network.lines.loc[[l[1] for l in line_objs]]
aggregated = {}
attr_inv = ["x", "r"]
attr_sum = ["s_nom", "b", "g", "s_nom_max", "s_nom_min"]
attr_mean = ["capital_cost", "length", "terrain_factor"]
for attr in attr_inv:
aggregated[attr] = 1.0 / (1.0 / lines[attr]).sum()
for attr in attr_sum:
aggregated[attr] = lines[attr].sum()
for attr in attr_mean:
aggregated[attr] = lines[attr].mean()
count += len(line_objs) - 1
# remove all but first line
for line in line_objs[1:]:
network.remove("Line", line[1])
rep = line_objs[0]
for key, value in aggregated.items():
setattr(rep, key, value)
seen.append((u, v))
logger.info(
"Removed %d excess lines from sub-network %s and replaced with aggregated lines",
count,
sub_network.name,
)
def find_tree(sub_network, weight="x_pu"):
"""
Get the spanning tree of the graph, choose the node with the highest degree
as a central "tree slack" and then see for each branch which paths from the
slack to each node go through the branch.
"""
branches_bus0 = sub_network.branches()["bus0"]
branches_i = branches_bus0.index
buses_i = sub_network.buses_i()
graph = sub_network.graph(weight=weight, inf_weight=1.0)
sub_network.tree = nx.minimum_spanning_tree(graph)
# find bus with highest degree to use as slack
tree_slack_bus, slack_degree = max(degree(sub_network.tree), key=itemgetter(1))
logger.debug("Tree slack bus is %s with degree %d.", tree_slack_bus, slack_degree)
# determine which buses are supplied in tree through branch from slack
# matrix to store tree structure
sub_network.T = dok_matrix((len(branches_i), len(buses_i)))
for j, bus in enumerate(buses_i):
path = nx.shortest_path(sub_network.tree, bus, tree_slack_bus)
for i in range(len(path) - 1):
branch = next(iter(graph[path[i]][path[i + 1]].keys()))
branch_i = branches_i.get_loc(branch)
sign = +1 if branches_bus0.iat[branch_i] == path[i] else -1
sub_network.T[branch_i, j] = sign
def find_cycles(sub_network, weight="x_pu"):
"""
Find all cycles in the sub_network and record them in sub_network.C.
networkx collects the cycles with more than 2 edges; then the 2-edge cycles
from the MultiGraph must be collected separately (for cases where there
are multiple lines between the same pairs of buses).
Cycles with infinite impedance are skipped.
"""
branches_bus0 = sub_network.branches()["bus0"]
branches_i = branches_bus0.index
# reduce to a non-multi-graph for cycles with > 2 edges
mgraph = sub_network.graph(weight=weight, inf_weight=False)
graph = nx.Graph(mgraph)
cycles = nx.cycle_basis(graph)
# number of 2-edge cycles
num_multi = len(mgraph.edges()) - len(graph.edges())
sub_network.C = dok_matrix((len(branches_bus0), len(cycles) + num_multi))
for j, cycle in enumerate(cycles):
for i in range(len(cycle)):
branch = next(iter(mgraph[cycle[i]][cycle[(i + 1) % len(cycle)]].keys()))
branch_i = branches_i.get_loc(branch)
sign = +1 if branches_bus0.iat[branch_i] == cycle[i] else -1
sub_network.C[branch_i, j] += sign
# counter for multis
c = len(cycles)
# add multi-graph 2-edge cycles for multiple branches between same pairs of buses
for u, v in graph.edges():
bs = list(mgraph[u][v].keys())
if len(bs) > 1:
first = bs[0]
first_i = branches_i.get_loc(first)
for b in bs[1:]:
b_i = branches_i.get_loc(b)
sign = (
-1 if branches_bus0.iat[b_i] == branches_bus0.iat[first_i] else +1
)
sub_network.C[first_i, c] = 1
sub_network.C[b_i, c] = sign
c += 1
def sub_network_lpf(sub_network, snapshots=None, skip_pre=False):
"""
Linear power flow for connected sub-network.
Parameters
----------
snapshots : list-like|single snapshot
A subset or an elements of network.snapshots on which to run
the power flow, defaults to network.snapshots
skip_pre : bool, default False
Skip the preliminary steps of computing topology, calculating
dependent values and finding bus controls.
Returns
-------
None
"""
snapshots = _as_snapshots(sub_network.network, snapshots)
logger.info(
"Performing linear load-flow on %s sub-network %s for snapshot(s) %s",
sub_network.network.sub_networks.at[sub_network.name, "carrier"],
sub_network,
snapshots,
)
network = sub_network.network
if not skip_pre:
calculate_dependent_values(network)
find_bus_controls(sub_network)
_allocate_pf_outputs(network, linear=True)
# get indices for the components on this subnetwork
buses_o = sub_network.buses_o
branches_i = sub_network.branches_i()
# allow all shunt impedances to dispatch as set
shunt_impedances_i = sub_network.shunt_impedances_i()
network.shunt_impedances_t.p.loc[
snapshots, shunt_impedances_i
] = network.shunt_impedances.g_pu.loc[shunt_impedances_i].values
# allow all one ports to dispatch as set
for c in sub_network.iterate_components(network.controllable_one_port_components):
c_p_set = get_switchable_as_dense(network, c.name, "p_set", snapshots, c.ind)
network.pnl(c.name).p.loc[snapshots, c.ind] = c_p_set
# set the power injection at each node
network.buses_t.p.loc[snapshots, buses_o] = sum(
[
(
(c.pnl.p.loc[snapshots, c.ind] * c.df.loc[c.ind, "sign"])
.groupby(c.df.loc[c.ind, "bus"], axis=1)
.sum()
.reindex(columns=buses_o, fill_value=0.0)
)
for c in sub_network.iterate_components(network.one_port_components)
]
+ [
(
-c.pnl["p" + str(i)]
.loc[snapshots]
.groupby(c.df["bus" + str(i)], axis=1)
.sum()
.reindex(columns=buses_o, fill_value=0)
)
for c in network.iterate_components(network.controllable_branch_components)
for i in [int(col[3:]) for col in c.df.columns if col[:3] == "bus"]
]
)
if not skip_pre and len(branches_i) > 0:
calculate_B_H(sub_network, skip_pre=True)
v_diff = np.zeros((len(snapshots), len(buses_o)))
if len(branches_i) > 0:
p = (
network.buses_t["p"].loc[snapshots, buses_o].values
- sub_network.p_bus_shift
)
v_diff[:, 1:] = spsolve(sub_network.B[1:, 1:], p[:, 1:].T).T
flows = (
pd.DataFrame(v_diff * sub_network.H.T, columns=branches_i, index=snapshots)
+ sub_network.p_branch_shift
)
for c in sub_network.iterate_components(network.passive_branch_components):
f = flows.loc[:, c.name]
network.pnl(c.name).p0.loc[snapshots, f.columns] = f
network.pnl(c.name).p1.loc[snapshots, f.columns] = -f
if network.sub_networks.at[sub_network.name, "carrier"] == "DC":
network.buses_t.v_mag_pu.loc[snapshots, buses_o] = 1 + v_diff
network.buses_t.v_ang.loc[snapshots, buses_o] = 0.0
else:
network.buses_t.v_ang.loc[snapshots, buses_o] = v_diff
network.buses_t.v_mag_pu.loc[snapshots, buses_o] = 1.0
# set slack bus power to pick up remained
slack_adjustment = (
-network.buses_t.p.loc[snapshots, buses_o[1:]].sum(axis=1).fillna(0.0)
- network.buses_t.p.loc[snapshots, buses_o[0]]
)
network.buses_t.p.loc[snapshots, buses_o[0]] += slack_adjustment
# let slack generator take up the slack
if sub_network.slack_generator is not None:
network.generators_t.p.loc[
snapshots, sub_network.slack_generator
] += slack_adjustment
def network_batch_lpf(network, snapshots=None):
"""
Batched linear power flow with numpy.dot for several snapshots.
"""
raise NotImplementedError("Batch linear power flow not supported yet.")
