import math
import numpy as np
import pandas as pd
import pypsa
from pypower.idx_brch import * # noqa: F403
from pypower.idx_bus import * # noqa: F403
from pypower.idx_cost import * # noqa: F403
from pypower.idx_gen import * # noqa: F403
[docs]
def to_powermodels(pypsa_net):
"""
Convert pypsa network to network dictionary format, using the pypower
structure as an intermediate steps
powermodels network dictionary:
https://lanl-ansi.github.io/PowerModels.jl/stable/network-data/
pypower caseformat:
https://github.com/rwl/PYPOWER/blob/master/pypower/caseformat.py
:param pypsa_net:
:return:
"""
# calculate per unit values
pypsa.pf.calculate_dependent_values(pypsa_net)
# convert pypsa network to pypower datastructure
ppc, loads_t, gens_t = pypsa2ppc(pypsa_net)
# conver pypower datastructure to powermodels network dictionary
pm = ppc2pm(ppc, pypsa_net)
return pm, loads_t, gens_t
[docs]
def convert_storage_series(timeseries):
if len(timeseries) == 0:
return {}
else:
storage = {"time_horizon": len(timeseries), "storage_data": {}}
for i, v in enumerate(timeseries.values):
storage["storage_data"][i + 1] = {"p_req": v}
return storage
# FIXME: Static storage data is exported from the eDisGo network rather than the PyPSA
# network as the capacity of network doesn't seem to be available there. For
# consistency with the rest of the conversion, it should be converted from PyPSA as
# well.
# TODO: This will (probably) not work if there are multiple storage units connected to
# the same bus.
[docs]
def add_storage_from_edisgo(edisgo_obj, psa_net, pm_dict):
"""
Read static storage data (position and capacity) from eDisGo and export to
Powermodels dict
"""
# Drop values that are not available
storage = pd.DataFrame(edisgo_obj.topology.storage_units_df[["bus", "p_nom"]])
# Rename storage parameters to PowerModels naming convention
storage.columns = ["storage_bus", "energy_rating"]
# Fill in missing values so that PowerModels doesn't complain
# TODO: Get (some of) these from eDisGo/PyPSA as well
storage["energy"] = 0.0
storage["charge_rating"] = 1.0
storage["discharge_rating"] = 1.0
storage["charge_efficiency"] = 0.9487
storage["discharge_efficiency"] = 0.9487
storage["ps"] = 0.0
storage["qs"] = 0.0
storage["qmin"] = 0.0
storage["qmax"] = 0.0
storage["r"] = 0.0
storage["x"] = 0.0
storage["p_loss"] = 0.0
storage["q_loss"] = 0.0
storage["standby_loss"] = 0.0
storage["status"] = 1
# Get Bus indices from PyPSA net
storage["storage_bus"] = [
psa_net.buses.index.get_loc(bus) + 1 for bus in storage["storage_bus"]
]
storage.index = [i + 1 for i in range(len(storage))]
# Add dedicated 'index' column because PowerModels likes it
storage["index"] = storage.index.to_list()
# Add to PowerModels statics Dict
pm_dict["storage"] = storage.to_dict(orient="index")
[docs]
def pypsa2ppc(psa_net):
"""Converter from pypsa data structure to pypower data structure
adapted from pandapower's pd2ppc converter
https://github.com/e2nIEE/pandapower/blob/911f300a96ee0ac062d82f7684083168ff052586/pandapower/pd2ppc.py
"""
# build static pypower structure
ppc = _init_ppc()
ppc["name"] = psa_net.name
_build_bus(psa_net, ppc)
_build_gen(psa_net, ppc)
_build_branch(psa_net, ppc)
_build_transformers(psa_net, ppc)
_build_load(psa_net, ppc)
# TODO STORAGE UNITS
_build_storage_units(psa_net, ppc)
# built dictionaries if timeseries is used
time_horizon = len(psa_net.loads_t["p_set"])
try:
load_dict = _build_load_dict(psa_net, ppc)
print(
"Dictionary for load timeseries of timehorizon {} created".format(
time_horizon
)
)
except IndexError:
print("No load timeseries. Create empty dicts for timeseries of load")
load_dict = dict()
try:
gen_dict = _build_generator_dict(psa_net, ppc)
print(
"Dictionary for generator timeseries of timehorizon {} created".format(
time_horizon
)
)
except IndexError:
print(
"no generator timeseries Create empty dicts "
"for timeseries of load and generation "
)
gen_dict = dict()
return ppc, load_dict, gen_dict
[docs]
def ppc2pm(ppc, psa_net): # pragma: no cover
"""
converter from pypower datastructure to powermodels dictionary,
adapted from pandapower to powermodels converter:
https://github.com/e2nIEE/pandapower/blob/develop/pandapower/converter/pandamodels/to_pm.py
:param ppc:
:return:
"""
pm = {
"gen": dict(),
"branch": dict(),
"bus": dict(),
"dcline": dict(),
"load": dict(),
"storage": dict(),
"baseMVA": ppc["baseMVA"],
"source_version": "2.0.0",
"shunt": dict(),
"sourcetype": "matpower",
"per_unit": True,
"name": ppc["name"],
}
load_idx = 1
shunt_idx = 1
for row in ppc["bus"]:
bus = dict()
idx = int(row[BUS_I]) + 1 # noqa: F405
bus["index"] = idx
bus["bus_i"] = idx
bus["zone"] = int(row[ZONE]) # noqa: F405
bus["bus_type"] = int(row[BUS_TYPE]) # noqa: F405
bus["vmax"] = row[VMAX] # noqa: F405
bus["vmin"] = row[VMIN] # noqa: F405
bus["va"] = row[VA] # noqa: F405
bus["vm"] = row[VM] # noqa: F405
bus["base_kv"] = row[BASE_KV] # noqa: F405
pd = row[PD] # noqa: F405
qd = row[QD] # noqa: F405
if pd != 0 or qd != 0:
pm["load"][str(load_idx)] = {
"pd": pd,
"qd": qd,
"load_bus": idx,
"status": True,
"index": load_idx,
}
load_idx += 1
bs = row[BS] # noqa: F405
gs = row[GS] # noqa: F405
if pd != 0 or qd != 0:
pm["shunt"][str(shunt_idx)] = {
"gs": gs,
"bs": bs,
"shunt_bus": idx,
"status": True,
"index": shunt_idx,
}
shunt_idx += 1
pm["bus"][str(idx)] = bus
n_lines = len(ppc["branch"])
for idx, row in enumerate(ppc["branch"], start=1):
branch = dict()
branch["index"] = idx
branch["transformer"] = idx > n_lines
branch["br_r"] = row[BR_R].real # noqa: F405
branch["br_x"] = row[BR_X].real # noqa: F405
branch["g_fr"] = -row[BR_B].imag / 2.0 # noqa: F405
branch["g_to"] = -row[BR_B].imag / 2.0 # noqa: F405
branch["b_fr"] = row[BR_B].real / 2.0 # noqa: F405
branch["b_to"] = row[BR_B].real / 2.0 # noqa: F405
branch["rate_a"] = (
row[RATE_A].real if row[RATE_A] > 0 else row[RATE_B].real # noqa: F405
)
branch["rate_b"] = row[RATE_B].real # noqa: F405
branch["rate_c"] = row[RATE_C].real # noqa: F405
branch["f_bus"] = int(row[F_BUS].real) + 1 # noqa: F405
branch["t_bus"] = int(row[T_BUS].real) + 1 # noqa: F405
branch["br_status"] = int(row[BR_STATUS].real) # noqa: F405
branch["angmin"] = row[ANGMIN].real # noqa: F405
branch["angmax"] = row[ANGMAX].real # noqa: F405
branch["tap"] = row[TAP].real # noqa: F405
branch["shift"] = math.radians(row[SHIFT].real) # noqa: F405
pm["branch"][str(idx)] = branch
for idx, row in enumerate(ppc["gen"], start=1):
gen = dict()
gen["pg"] = row[PG] # noqa: F405
gen["qg"] = row[QG] # noqa: F405
gen["gen_bus"] = int(row[GEN_BUS]) + 1 # noqa: F405
gen["vg"] = row[VG] # noqa: F405
gen["qmax"] = row[QMAX] # noqa: F405
gen["gen_status"] = int(row[GEN_STATUS]) # noqa: F405
gen["qmin"] = row[QMIN] # noqa: F405
gen["pmin"] = row[PMIN] # noqa: F405
gen["pmax"] = row[PMAX] # noqa: F405
gen["index"] = idx
pm["gen"][str(idx)] = gen
# TODO add attribute "fluctuating" to generators from psa_net, maybe move to ppc
# first
# is_fluctuating = [
# int("fluctuating" in index.lower()) for index in psa_net.generators.index
# ]
# for idx, row in enumerate(is_fluctuating, start=1):
# pm["gen"][str(idx)]["fluctuating"] = row
for idx, row in enumerate(psa_net.generators["control"], start=1):
pm["gen"][str(idx)]["gen_slack"] = (row == "Slack") * 1
for idx, row in enumerate(psa_net.generators["fluctuating"], start=1):
# convert boolean to 0 and 1, check if row is nan, e.g. slack bus
pm["gen"][str(idx)]["fluctuating"] = (not (math.isnan(row)) and row) * 1
if len(ppc["gencost"]) > len(ppc["gen"]):
ppc["gencost"] = ppc["gencost"][: ppc["gen"].shape[0], :]
for idx, row in enumerate(ppc["gencost"], start=1):
gen = pm["gen"][str(idx)]
gen["model"] = int(row[MODEL]) # noqa: F405
if gen["model"] == 1:
gen["ncost"] = int(row[NCOST]) # noqa: F405
gen["cost"] = row[COST : COST + gen["ncost"] * 2].tolist() # noqa: F405
elif gen["model"] == 2:
gen["ncost"] = int(row[NCOST]) # noqa: F405
gen["cost"] = [0] * 3
costs = row[COST:] # noqa: F405
if len(costs) > 3:
print(costs)
raise ValueError("Maximum quadratic cost function allowed")
gen["cost"][-len(costs) :] = costs
if len(ppc["branchcost"]) > len(ppc["branch"]):
ppc["branchcost"] = ppc["branchcost"][: ppc["branch"].shape[0], :]
for idx, row in enumerate(ppc["branchcost"], start=1):
ncost = int(row[0])
branch = pm["branch"][str(idx)]
branch["ncost"] = ncost
branch["cost"] = [0] * ncost
costs = row[1:]
branch["cost"][-len(costs) :] = costs
# TODO STORAGE UNITS!
return pm
def _init_ppc():
# init empty ppc
ppc = {
"name": "name",
"baseMVA": 1,
"version": 2,
"bus": np.array([], dtype=float),
"branch": np.array([], dtype=np.complex128),
"gen": np.array([], dtype=float),
"gencost": np.array([], dtype=float),
"branchcost": np.array([], dtype=float),
"internal": {
"Ybus": np.array([], dtype=np.complex128),
"Yf": np.array([], dtype=np.complex128),
"Yt": np.array([], dtype=np.complex128),
"branch_is": np.array([], dtype=bool),
"gen_is": np.array([], dtype=bool),
"DLF": np.array([], dtype=np.complex128),
"buses_ord_bfs_nets": np.array([], dtype=float),
},
}
return ppc
def _build_bus(psa_net, ppc):
n_bus = len(psa_net.buses.index)
print("build {} buses".format(n_bus))
col_names = (
"index",
"type",
"Pd",
"Qd",
"Gs",
"Bs",
"area",
"v_mag_pu_set",
"v_ang_set",
"v_nom",
"zone",
"v_mag_pu_max",
"v_mag_pu_min".split(", "),
)
bus_cols = len(col_names)
ppc["bus"] = np.zeros(shape=(n_bus, bus_cols), dtype=float)
ppc["bus"][:, :bus_cols] = np.array([0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1.05, 0.95])
ppc["bus"][:, BUS_I] = np.arange(n_bus) # noqa: F405
bus_types = ["PQ", "PV", "Slack", "None"]
bus_types_int = np.array(
[bus_types.index(b_type) + 1 for b_type in psa_net.buses["control"].values],
dtype=int,
)
ppc["bus"][:, BUS_TYPE] = bus_types_int # noqa: F405
ppc["bus"][:, BASE_KV] = psa_net.buses["v_nom"].values # noqa: F405
# for edisgo scenario voltage bounds defined for load and feed-in case with
# 0.985<= v <= 1.05 bounds have to be at least in that range, only accept stronger
# bounds if given
ppc["bus"][:, VMAX] = [ # noqa: F405
min(val, 1.05) for val in psa_net.buses["v_mag_pu_max"].values
]
ppc["bus"][:, VMIN] = [ # noqa: F405
max(val, 0.985) for val in psa_net.buses["v_mag_pu_min"].values
]
return
def _build_gen(psa_net, ppc):
n_gen = psa_net.generators.shape[0]
gen_cols = 21
# "bus, p_set, q_set, q_max, q_min, v_set_pu, mva_base, status, p_nom, p_min, Pc1,
# Pc2, Qc1min, Qc1max, Qc2min, Qc2max, ramp_agc, ramp_10, ramp_30, ramp_q, apf
ppc["gen"] = np.zeros(shape=(n_gen, gen_cols), dtype=float)
# get bus indices for generators
bus_indices = np.array(
[
psa_net.buses.index.get_loc(bus_name)
for bus_name in psa_net.generators["bus"]
]
)
print(
"build {} generators, distributed on {} buses".format(
n_gen, len(np.unique(bus_indices))
)
)
ppc["gen"][:, GEN_BUS] = bus_indices # noqa: F405
# adjust bus types
# bus_types = ["PQ", "PV", "Slack", "None"]
# gen_types = np.array(
# [
# bus_types.index(gen_type) + 1
# for gen_type in psa_net.generators["control"].values
# ],
# dtype=int,
# )
# ppc["bus"][bus_indices,BUS_TYPE] = gen_types
# set setpoint of pg and qg
ppc["gen"][:, PG] = psa_net.generators["p_set"].values # noqa: F405
ppc["gen"][:, QG] = psa_net.generators["q_set"].values # noqa: F405
ppc["gen"][:, MBASE] = 1.0 # noqa: F405
ppc["gen"][:, GEN_STATUS] = 1.0 # noqa: F405
ppc["gen"][:, PMAX] = psa_net.generators["p_nom"].values # noqa: F405
ppc["gen"][:, PMIN] = 0 # noqa: F405
# TODO SET QMAX AND QMIN! e.g.: cos(phi) value from config
# ppc["gen"][:,QMAX] = 0
# ppc["gen"][:, QMIN] = 0
# build field for generator costs
# 2 startup shutdown n c(n-1) ... c0
# for quadratic cost function n=3--> c2,c1,c0 lead to 7 columns
cost_cols = 7
ppc["gencost"] = np.zeros(shape=(n_gen, cost_cols), dtype=float)
# polynomial cost function
ppc["gencost"][:, MODEL] = POLYNOMIAL # noqa: F405
ppc["gencost"][:, STARTUP] = psa_net.generators[ # noqa: F405
"start_up_cost"
].values
ppc["gencost"][:, SHUTDOWN] = psa_net.generators[ # noqa: F405
"shut_down_cost"
].values
# quadratic cost function has 3 cost coefficients
ppc["gencost"][:, NCOST] = 3 # noqa: F405
ppc["gencost"][:, COST] = 0.0 # noqa: F405
ppc["gencost"][:, COST + 1] = psa_net.generators[ # noqa: F405
"marginal_cost"
].values
ppc["gencost"][:, COST + 2] = 0.0 # noqa: F405
return
def _build_branch(psa_net, ppc):
n_branch = len(psa_net.lines.index)
print("build {} lines".format(n_branch))
col_names = [
"fbus",
"tbus",
"r",
"x",
"b",
"rateA",
"rateB",
"rateC",
"ratio",
"angle",
"status",
"angmin",
"angmax",
]
branch_cols = len(col_names)
ppc["branch"] = np.zeros(shape=(n_branch, branch_cols), dtype=float)
from_bus = np.array(
[psa_net.buses.index.get_loc(bus_name) for bus_name in psa_net.lines["bus0"]]
)
to_bus = np.array(
[psa_net.buses.index.get_loc(bus_name) for bus_name in psa_net.lines["bus1"]]
)
ppc["branch"][:, F_BUS] = from_bus # noqa: F405
ppc["branch"][:, T_BUS] = to_bus # noqa: F405
ppc["branch"][:, BR_R] = psa_net.lines["r_pu"].values # noqa: F405
ppc["branch"][:, BR_X] = psa_net.lines["x_pu"].values # noqa: F405
ppc["branch"][:, BR_B] = psa_net.lines["b_pu"].values # noqa: F405
ppc["branch"][:, RATE_A] = psa_net.lines["s_nom"].values # noqa: F405
ppc["branch"][:, RATE_B] = 250 # Default values # noqa: F405
ppc["branch"][:, RATE_C] = 250 # Default values # noqa: F405
ppc["branch"][:, TAP] = 0.0 # noqa: F405
ppc["branch"][:, SHIFT] = 0.0 # noqa: F405
ppc["branch"][:, BR_STATUS] = 1.0 # noqa: F405
ppc["branch"][:, ANGMIN] = -360 # noqa: F405
ppc["branch"][:, ANGMAX] = 360 # noqa: F405
# TODO BRANCHCOSTS!
# check which branch costs are given in psa_net,
ncost = sum(
(colName in psa_net.lines.columns) * 1
for colName in ["costs_earthworks", "costs_cable"]
)
if ncost == 0:
print("no branch costs are given in pypsa network")
elif ncost == 1:
if "costs_cable" not in psa_net.lines.columns:
print(
"costs for cables not in pypsa network, not possible to define cost"
"function for network expansion"
)
else:
ppc["branchcost"] = np.zeros(shape=(n_branch, 2), dtype=float)
ppc["branchcost"][:, 0] = ncost
ppc["branchcost"][:, 1] = psa_net.lines["costs_cable"].values
elif ncost == 2:
ppc["branchcost"] = np.zeros(shape=(n_branch, 3), dtype=float)
ppc["branchcost"][:, 0] = ncost
ppc["branchcost"][:, 1] = psa_net.lines["costs_cable"].values
ppc["branchcost"][:, 2] = psa_net.lines["costs_earthworks"].values
return
def _build_transformers(psa_net, ppc):
n_transformers = len(psa_net.transformers.index)
print("appending {} transformers".format(n_transformers))
col_names = [
"fbus",
"tbus",
"r",
"x",
"b",
"rateA",
"rateB",
"rateC",
"ratio",
"angle",
"status",
"angmin",
"angmax",
]
transformers = np.zeros(shape=(n_transformers, len(col_names)), dtype=float)
from_bus = np.array(
[
psa_net.buses.index.get_loc(bus_name)
for bus_name in psa_net.transformers["bus0"]
]
)
to_bus = np.array(
[
psa_net.buses.index.get_loc(bus_name)
for bus_name in psa_net.transformers["bus1"]
]
)
transformers[:, F_BUS] = from_bus # noqa: F405
transformers[:, T_BUS] = to_bus # noqa: F405
transformers[:, BR_R] = psa_net.transformers["r_pu"].values # noqa: F405
transformers[:, BR_X] = psa_net.transformers["x_pu"].values # noqa: F405
transformers[:, BR_B] = psa_net.transformers["b_pu"].values # noqa: F405
transformers[:, RATE_A] = psa_net.transformers["s_nom"].values # noqa: F405
transformers[:, RATE_B] = 250 # Default values # noqa: F405
transformers[:, RATE_C] = 250 # Default values # noqa: F405
transformers[:, TAP] = psa_net.transformers["tap_ratio"].values # noqa: F405
transformers[:, SHIFT] = psa_net.transformers["phase_shift"].values # noqa: F405
transformers[:, BR_STATUS] = 1.0 # noqa: F405
transformers[:, ANGMIN] = -360 # noqa: F405
transformers[:, ANGMAX] = 360 # noqa: F405
ppc["branch"] = np.append(ppc["branch"], transformers, axis=0)
# add trafo costs to branch cost with same shape
if len(ppc["branchcost"]) > 0:
print("append transformer costs")
ncost = ppc["branchcost"].shape[1] - 1
trafo_costs = np.zeros(shape=(n_transformers, ncost + 1), dtype=float)
if hasattr(psa_net.transformers, "trafo_costs"):
trafo_costs[:, 0] = ncost
trafo_costs[:, 1] = psa_net.transformers["trafo_costs"].values
print(trafo_costs)
ppc["branchcost"] = np.append(ppc["branchcost"], trafo_costs, axis=0)
return
def _build_load(psa_net, ppc):
n_load = psa_net.loads.shape[0]
load_buses = np.array(
[psa_net.buses.index.get_loc(bus_name) for bus_name in psa_net.loads["bus"]]
)
print(
"build {} loads, distributed on {} buses".format(
n_load, len(np.unique(load_buses))
)
)
# USE LOAD DATA FROM psa_net.loads as static network data
# set bool if loads contains a timeseries
# istime = len(psa_net.loads_t["p_set"].values[0]) != 0
# istime = False
# print("network has timeseries for load: {}".format(istime))
for load_idx, bus_idx in enumerate(load_buses):
# if istime:
# # if timeseries take maximal value of load_bus for static information of
# # the network
# p_d = max(psa_net.loads_t["p_set"].values[:,load_idx])
# q_d = max(psa_net.loads_t["q_set"].values[:,load_idx])
# else:
p_d = psa_net.loads["p_set"].values[load_idx]
q_d = psa_net.loads["q_set"].values[load_idx]
# increase demand at bus_idx by p_d and q_d from load_idx, as multiple loads
# can be attached to single bus
ppc["bus"][bus_idx, PD] += p_d # noqa: F405
ppc["bus"][bus_idx, QD] += q_d # noqa: F405
return
def _build_storage_units(psa_net, ppc):
print("storage units are not implemented yet")
def _build_load_dict(psa_net, ppc):
"""
build load dict containing timeseries from psa_net.loads_t
:param psa_net: pypsa network
:param ppc:
:return: load_dict: Dict()
"""
load_dict = {"load_data": dict()}
load_buses = np.array(
[psa_net.buses.index.get_loc(bus_name) for bus_name in psa_net.loads["bus"]]
)
time_horizon = len(psa_net.loads_t["p_set"])
load_dict["time_horizon"] = time_horizon
for t in range(time_horizon):
load_dict["load_data"][str(t + 1)] = dict()
for load_idx, bus_idx in enumerate(load_buses):
# p_d = psa_net.loads_t["p_set"].values[t,load_idx]
# qd = psa_net.loads_t["q_set"].values[t,load_idx]
p_d = psa_net.loads_t["p_set"][psa_net.loads.index[load_idx]][t]
qd = psa_net.loads_t["q_set"][psa_net.loads.index[load_idx]][t]
load_dict["load_data"][str(t + 1)][str(load_idx + 1)] = {
"pd": p_d,
"qd": qd,
"load_bus": int(bus_idx + 1),
"status": True,
"index": int(load_idx + 1),
}
return load_dict
def _build_generator_dict(psa_net, ppc):
generator_dict = {"gen_data": dict()}
time_horizon = len(psa_net.generators_t["p_set"])
generator_dict["time_horizon"] = time_horizon
# buses_with_gens = [psa_net.generators.loc[busname]["bus"] for busname in
# psa_net.generators_t["p_set"].columns]
# gen_buses = np.array([psa_net.buses.index.get_loc(bus_name) for bus_name in
# buses_with_gens])
gen_buses = [
psa_net.buses.index.get_loc(bus_name) for bus_name in psa_net.generators["bus"]
]
for t in range(time_horizon):
generator_dict["gen_data"][str(t + 1)] = dict()
for gen_idx, bus_idx in enumerate(gen_buses):
# pg = psa_net.generators_t["p_set"].values[t, gen_idx]
# qg = psa_net.generators_t["q_set"].values[t, gen_idx]
pg = psa_net.generators_t["p_set"][psa_net.generators.index[gen_idx]][t]
qg = psa_net.generators_t["q_set"][psa_net.generators.index[gen_idx]][t]
# if no value is set, set pg and qg to large value, e.g. representing slack
# TODO verify or find another solution not using "large" value
if np.isnan(pg):
pg = 99999
if np.isnan(qg):
qg = 99999
generator_dict["gen_data"][str(t + 1)][str(gen_idx + 1)] = {
"pg": pg,
"qg": qg,
"gen_bus": int(bus_idx) + 1,
"status": True,
"index": int(gen_idx + 1),
}
return generator_dict