e2nIEE / pandapower
 1 ```# -*- coding: utf-8 -*- ``` 2 3 ```# Copyright (c) 2016-2021 by University of Kassel and Fraunhofer Institute for Energy Economics ``` 4 ```# and Energy System Technology (IEE), Kassel. All rights reserved. ``` 5 6 7 1 ```from collections import defaultdict, Counter ``` 8 1 ```from itertools import chain ``` 9 10 1 ```import numpy as np ``` 11 1 ```import pandas as pd ``` 12 13 1 ```from pandapower.auxiliary import _sum_by_group ``` 14 1 ```from pandapower.pypower.idx_bus import BUS_I, BASE_KV, PD, QD, GS, BS, VMAX, VMIN, BUS_TYPE, NONE, VM, VA, \ ``` 15 ``` CID, CZD, bus_cols, REF ``` 16 17 1 ```try: ``` 18 1 ``` from numba import jit ``` 19 0 ```except ImportError: ``` 20 0 ``` from .pf.no_numba import jit ``` 21 22 23 1 ```@jit(nopython=True, cache=False) ``` 24 ```def ds_find(ar, bus): # pragma: no cover ``` 25 ``` while True: ``` 26 ``` p = ar[bus] ``` 27 ``` if p == bus: ``` 28 ``` break ``` 29 ``` bus = p ``` 30 ``` return p ``` 31 32 33 1 ```@jit(nopython=True, cache=False) ``` 34 ```def ds_union(ar, bus1, bus2, bus_is_pv): # pragma: no cover ``` 35 ``` root1 = ds_find(ar, bus1) ``` 36 ``` root2 = ds_find(ar, bus2) ``` 37 ``` if root1 == root2: ``` 38 ``` return ``` 39 ``` if bus_is_pv[root2]: ``` 40 ``` ar[root1] = root2 ``` 41 ``` else: ``` 42 ``` ar[root2] = root1 ``` 43 44 45 1 ```@jit(nopython=True, cache=False) ``` 46 ```def ds_create(ar, switch_bus, switch_elm, switch_et_bus, switch_closed, switch_z_ohm, ``` 47 ``` bus_is_pv, bus_in_service): # pragma: no cover ``` 48 ``` for i in range(len(switch_bus)): ``` 49 ``` if not switch_closed[i] or not switch_et_bus[i] or switch_z_ohm[i] > 0: ``` 50 ``` continue ``` 51 ``` bus1 = switch_bus[i] ``` 52 ``` bus2 = switch_elm[i] ``` 53 ``` if bus_in_service[bus1] and bus_in_service[bus2]: ``` 54 ``` ds_union(ar, bus1, bus2, bus_is_pv) ``` 55 56 57 1 ```@jit(nopython=True, cache=False) ``` 58 1 ```def fill_bus_lookup(ar, bus_lookup, bus_index): ``` 59 0 ``` for i in range(len(bus_index)): ``` 60 0 ``` bus_lookup[bus_index[i]] = i ``` 61 0 ``` for b in bus_index: ``` 62 0 ``` ds = ds_find(ar, b) ``` 63 0 ``` bus_lookup[b] = bus_lookup[ar[ds]] ``` 64 65 66 1 ```def create_bus_lookup_numba(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask): ``` 67 1 ``` max_bus_idx = np.max(bus_index) ``` 68 ``` # extract numpy arrays of switch table data ``` 69 1 ``` switch = net["switch"] ``` 70 1 ``` switch_bus = switch["bus"].values ``` 71 1 ``` switch_elm = switch["element"].values ``` 72 1 ``` switch_et_bus = switch["et"].values == "b" ``` 73 1 ``` switch_closed = switch["closed"].values ``` 74 1 ``` switch_z_ohm = switch['z_ohm'].values ``` 75 ``` # create array for fast checking if a bus is in_service ``` 76 1 ``` bus_in_service = np.zeros(max_bus_idx + 1, dtype=bool) ``` 77 1 ``` bus_in_service[bus_is_idx] = True ``` 78 ``` # create array for fast checking if a bus is pv bus ``` 79 1 ``` bus_is_pv = np.zeros(max_bus_idx + 1, dtype=bool) ``` 80 1 ``` bus_is_pv[net["ext_grid"]["bus"].values[eg_is_mask]] = True ``` 81 1 ``` bus_is_pv[net["gen"]["bus"].values[gen_is_mask]] = True ``` 82 ``` # create array that represents the disjoint set ``` 83 1 ``` ar = np.arange(max_bus_idx + 1) ``` 84 1 ``` ds_create(ar, switch_bus, switch_elm, switch_et_bus, switch_closed, switch_z_ohm, bus_is_pv, bus_in_service) ``` 85 ``` # finally create and fill bus lookup ``` 86 1 ``` bus_lookup = -np.ones(max_bus_idx + 1, dtype=int) ``` 87 1 ``` fill_bus_lookup(ar, bus_lookup, bus_index) ``` 88 1 ``` return bus_lookup ``` 89 90 91 1 ```class DisjointSet(dict): ``` 92 1 ``` def add(self, item): ``` 93 0 ``` self[item] = item ``` 94 95 1 ``` def find(self, item): ``` 96 1 ``` parent = self[item] ``` 97 1 ``` if self[parent] != parent: ``` 98 1 ``` parent = self.find(parent) ``` 99 1 ``` self[item] = parent ``` 100 1 ``` return parent ``` 101 102 1 ``` def union(self, item1, item2): ``` 103 1 ``` p1 = self.find(item1) ``` 104 1 ``` p2 = self.find(item2) ``` 105 1 ``` self[p1] = p2 ``` 106 107 108 1 ```def create_consecutive_bus_lookup(net, bus_index): ``` 109 ``` # create a mapping from arbitrary pp-index to a consecutive index starting at zero (ppc-index) ``` 110 1 ``` consec_buses = np.arange(len(bus_index)) ``` 111 ``` # bus_lookup as dict: ``` 112 ``` # bus_lookup = dict(zip(bus_index, consec_buses)) ``` 113 ``` # bus lookup as mask from pandapower -> pypower ``` 114 1 ``` bus_lookup = -np.ones(max(bus_index) + 1, dtype=int) ``` 115 1 ``` bus_lookup[bus_index] = consec_buses ``` 116 1 ``` return bus_lookup ``` 117 118 119 1 ```def create_bus_lookup_numpy(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask, closed_bb_switch_mask): ``` 120 1 ``` bus_lookup = create_consecutive_bus_lookup(net, bus_index) ``` 121 1 ``` net._fused_bb_switches = closed_bb_switch_mask & (net["switch"]["z_ohm"].values <= 0) ``` 122 1 ``` if net._fused_bb_switches.any(): ``` 123 ``` # Note: this might seem a little odd - first constructing a pp to ppc mapping without ``` 124 ``` # fused busses and then update the entries. The alternative (to construct the final ``` 125 ``` # mapping at once) would require to determine how to fuse busses and which busses ``` 126 ``` # are not part of any opened bus-bus switches first. It turns out, that the latter takes ``` 127 ``` # quite some time in the average usecase, where #busses >> #bus-bus switches. ``` 128 129 ``` # Find PV / Slack nodes -> their bus must be kept when fused with a PQ node ``` 130 1 ``` pv_list = [net["ext_grid"]["bus"].values[eg_is_mask], net["gen"]["bus"].values[gen_is_mask]] ``` 131 1 ``` pv_ref = np.unique(np.hstack(pv_list)) ``` 132 ``` # get the pp-indices of the buses which are connected to a switch to be fused ``` 133 1 ``` fbus = net["switch"]["bus"].values[net._fused_bb_switches] ``` 134 1 ``` tbus = net["switch"]["element"].values[net._fused_bb_switches] ``` 135 136 ``` # create a mapping to map each bus to itself at frist ... ``` 137 1 ``` ds = DisjointSet({e: e for e in chain(fbus, tbus)}) ``` 138 139 ``` # ... to follow each bus along a possible chain of switches to its final bus and update the ``` 140 ``` # map ``` 141 1 ``` for f, t in zip(fbus, tbus): ``` 142 1 ``` ds.union(f, t) ``` 143 144 ``` # now we can find out how to fuse each bus by looking up the final bus of the chain they ``` 145 ``` # are connected to ``` 146 1 ``` v = defaultdict(set) ``` 147 1 ``` for a in ds: ``` 148 1 ``` v[ds.find(a)].add(a) ``` 149 ``` # build sets of buses which will be fused ``` 150 1 ``` disjoint_sets = [e for e in v.values() if len(e) > 1] ``` 151 152 ``` # check if PV buses need to be fused ``` 153 ``` # if yes: the sets with PV buses must be found (which is slow) ``` 154 ``` # if no: the check can be omitted ``` 155 1 ``` if any(i in fbus or i in tbus for i in pv_ref): ``` 156 ``` # for every disjoint set ``` 157 1 ``` for dj in disjoint_sets: ``` 158 ``` # check if pv buses are in the disjoint set dj ``` 159 1 ``` pv_buses_in_set = set(pv_ref) & dj ``` 160 1 ``` nr_pv_bus = len(pv_buses_in_set) ``` 161 1 ``` if nr_pv_bus == 0: ``` 162 ``` # no pv buses. Use any bus in dj ``` 163 1 ``` map_to = bus_lookup[dj.pop()] ``` 164 ``` else: ``` 165 ``` # one pv bus. Get bus from pv_buses_in_set ``` 166 1 ``` map_to = bus_lookup[pv_buses_in_set.pop()] ``` 167 1 ``` for bus in dj: ``` 168 ``` # update lookup ``` 169 1 ``` bus_lookup[bus] = map_to ``` 170 ``` else: ``` 171 ``` # no PV buses in set ``` 172 1 ``` for dj in disjoint_sets: ``` 173 ``` # use any bus in set ``` 174 1 ``` map_to = bus_lookup[dj.pop()] ``` 175 1 ``` for bus in dj: ``` 176 ``` # update bus lookup ``` 177 1 ``` bus_lookup[bus] = map_to ``` 178 1 ``` return bus_lookup ``` 179 180 181 1 ```def create_bus_lookup(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask, numba): ``` 182 1 ``` switches_with_pos_z_ohm = net["switch"]["z_ohm"].values > 0 ``` 183 1 ``` if switches_with_pos_z_ohm.any() or numba == False: ``` 184 ``` # if there are any closed bus-bus switches find them ``` 185 1 ``` closed_bb_switch_mask = ((net["switch"]["closed"].values == True) & ``` 186 ``` (net["switch"]["et"].values == "b") & ``` 187 ``` np.in1d(net["switch"]["bus"].values, bus_is_idx) & ``` 188 ``` np.in1d(net["switch"]["element"].values, bus_is_idx)) ``` 189 190 1 ``` if switches_with_pos_z_ohm.any(): ``` 191 1 ``` net._impedance_bb_switches = closed_bb_switch_mask & switches_with_pos_z_ohm ``` 192 ``` else: ``` 193 1 ``` net._impedance_bb_switches = np.zeros(switches_with_pos_z_ohm.shape) ``` 194 195 1 ``` if numba: ``` 196 1 ``` bus_lookup = create_bus_lookup_numba(net, bus_index, bus_is_idx, ``` 197 ``` gen_is_mask, eg_is_mask) ``` 198 ``` else: ``` 199 1 ``` bus_lookup = create_bus_lookup_numpy(net, bus_index, bus_is_idx, ``` 200 ``` gen_is_mask, eg_is_mask, closed_bb_switch_mask) ``` 201 1 ``` return bus_lookup ``` 202 203 204 1 ```def get_voltage_init_vector(net, init_v, mode): ``` 205 1 ``` if isinstance(init_v, str): ``` 206 1 ``` if init_v == "results": ``` 207 ``` # init voltage possible if bus results are available ``` 208 1 ``` if "res_bus" in net and net.res_bus.index.equals(net.bus.index): ``` 209 ``` # init bus voltages from results if the sorting is correct ``` 210 1 ``` res_table = "res_bus" ``` 211 ``` else: ``` 212 ``` # cannot init from results, since sorting of results is different from element table ``` 213 ``` # TO BE REVIEWED! Why there was no raise before this commit? ``` 214 0 ``` raise UserWarning("Init from results not possible. Index of res_bus do not match with bus. " ``` 215 ``` "You should sort res_bus before calling runpp.") ``` 216 0 ``` return None ``` 217 218 1 ``` if mode == "magnitude": ``` 219 1 ``` return net[res_table]["vm_pu"].values.copy() ``` 220 1 ``` elif mode == "angle": ``` 221 1 ``` return net[res_table]["va_degree"].values.copy() ``` 222 ``` else: ``` 223 0 ``` raise UserWarning(str(mode)+" for initialization not available!") ``` 224 1 ``` if init_v == "flat": ``` 225 1 ``` if mode == "magnitude": ``` 226 1 ``` net["_options"]["init_vm_pu"] = 0. ``` 227 ``` else: ``` 228 1 ``` net["_options"]["init_va_degree"] = 0. ``` 229 1 ``` return None ``` 230 1 ``` elif isinstance(init_v, (float, np.ndarray, list)): ``` 231 1 ``` return init_v ``` 232 1 ``` elif isinstance(init_v, pd.Series): ``` 233 0 ``` if init_v.index.equals(net.bus.index): ``` 234 0 ``` return init_v.loc[net.bus.index] ``` 235 ``` else: ``` 236 0 ``` raise UserWarning("Voltage starting vector indices do not match bus indices") ``` 237 238 239 1 ```def _build_bus_ppc(net, ppc): ``` 240 ``` """ ``` 241 ``` Generates the ppc["bus"] array and the lookup pandapower indices -> ppc indices ``` 242 ``` """ ``` 243 1 ``` init_vm_pu = net["_options"]["init_vm_pu"] ``` 244 1 ``` init_va_degree = net["_options"]["init_va_degree"] ``` 245 1 ``` mode = net["_options"]["mode"] ``` 246 1 ``` numba = net["_options"]["numba"] if "numba" in net["_options"] else False ``` 247 248 ``` # get bus indices ``` 249 1 ``` nr_xward = len(net.xward) ``` 250 1 ``` nr_trafo3w = len(net.trafo3w) ``` 251 1 ``` aux = dict() ``` 252 1 ``` if nr_xward > 0 or nr_trafo3w > 0: ``` 253 1 ``` bus_indices = [net["bus"].index.values, np.array([], dtype=np.int64)] ``` 254 1 ``` max_idx = max(net["bus"].index) + 1 ``` 255 1 ``` if nr_xward > 0: ``` 256 1 ``` aux_xward = np.arange(max_idx, max_idx + nr_xward, dtype=np.int64) ``` 257 1 ``` aux["xward"] = aux_xward ``` 258 1 ``` bus_indices.append(aux_xward) ``` 259 1 ``` if nr_trafo3w: ``` 260 1 ``` aux_trafo3w = np.arange(max_idx + nr_xward, max_idx + nr_xward + nr_trafo3w) ``` 261 1 ``` aux["trafo3w"] = aux_trafo3w ``` 262 1 ``` bus_indices.append(aux_trafo3w) ``` 263 1 ``` bus_index = np.concatenate(bus_indices) ``` 264 ``` else: ``` 265 1 ``` bus_index = net["bus"].index.values ``` 266 ``` # get in service elements ``` 267 268 1 ``` if mode == "nx": ``` 269 1 ``` bus_lookup = create_consecutive_bus_lookup(net, bus_index) ``` 270 ``` else: ``` 271 1 ``` _is_elements = net["_is_elements"] ``` 272 1 ``` eg_is_mask = _is_elements['ext_grid'] ``` 273 1 ``` gen_is_mask = _is_elements['gen'] ``` 274 1 ``` bus_is_idx = _is_elements['bus_is_idx'] ``` 275 1 ``` bus_lookup = create_bus_lookup(net, bus_index, bus_is_idx, ``` 276 ``` gen_is_mask, eg_is_mask, numba=numba) ``` 277 278 1 ``` n_bus_ppc = len(bus_index) ``` 279 ``` # init ppc with empty values ``` 280 1 ``` ppc["bus"] = np.zeros(shape=(n_bus_ppc, bus_cols), dtype=float) ``` 281 1 ``` ppc["bus"][:, :15] = np.array([0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 2, 0, 0., 0.]) # changes of ``` 282 ``` # voltage limits (2 and 0) must be considered in check_opf_data ``` 283 1 ``` if mode == "sc": ``` 284 1 ``` from pandapower.shortcircuit.idx_bus import bus_cols_sc ``` 285 1 ``` bus_sc = np.empty(shape=(n_bus_ppc, bus_cols_sc), dtype=float) ``` 286 1 ``` bus_sc.fill(np.nan) ``` 287 1 ``` ppc["bus"] = np.hstack((ppc["bus"], bus_sc)) ``` 288 289 ``` # apply consecutive bus numbers ``` 290 1 ``` ppc["bus"][:, BUS_I] = np.arange(n_bus_ppc) ``` 291 292 1 ``` n_bus = len(net.bus.index) ``` 293 ``` # init voltages from net ``` 294 1 ``` ppc["bus"][:n_bus, BASE_KV] = net["bus"]["vn_kv"].values ``` 295 ``` # set buses out of service (BUS_TYPE == 4) ``` 296 1 ``` if nr_xward > 0 or nr_trafo3w > 0: ``` 297 1 ``` in_service = np.concatenate([net["bus"]["in_service"].values, ``` 298 ``` net["xward"]["in_service"].values, ``` 299 ``` net["trafo3w"]["in_service"].values]) ``` 300 ``` else: ``` 301 1 ``` in_service = net["bus"]["in_service"].values ``` 302 1 ``` ppc["bus"][~in_service, BUS_TYPE] = NONE ``` 303 1 ``` if mode != "nx": ``` 304 1 ``` set_reference_buses(net, ppc, bus_lookup, mode) ``` 305 1 ``` vm_pu = get_voltage_init_vector(net, init_vm_pu, "magnitude") ``` 306 1 ``` if vm_pu is not None: ``` 307 1 ``` ppc["bus"][:n_bus, VM] = vm_pu ``` 308 309 1 ``` va_degree = get_voltage_init_vector(net, init_va_degree, "angle") ``` 310 1 ``` if va_degree is not None: ``` 311 1 ``` ppc["bus"][:n_bus, VA] = va_degree ``` 312 313 1 ``` if mode == "sc": ``` 314 1 ``` _add_c_to_ppc(net, ppc) ``` 315 316 1 ``` if net._options["mode"] == "opf": ``` 317 1 ``` if "max_vm_pu" in net.bus: ``` 318 1 ``` ppc["bus"][:n_bus, VMAX] = net["bus"].max_vm_pu.values ``` 319 ``` else: ``` 320 1 ``` ppc["bus"][:n_bus, VMAX] = 2. # changes of VMAX must be considered in check_opf_data ``` 321 1 ``` if "min_vm_pu" in net.bus: ``` 322 1 ``` ppc["bus"][:n_bus, VMIN] = net["bus"].min_vm_pu.values ``` 323 ``` else: ``` 324 1 ``` ppc["bus"][:n_bus, VMIN] = 0. # changes of VMIN must be considered in check_opf_data ``` 325 326 1 ``` if len(net.xward): ``` 327 1 ``` _fill_auxiliary_buses(net, ppc, bus_lookup, "xward", "bus", aux) ``` 328 329 1 ``` if len(net.trafo3w): ``` 330 1 ``` _fill_auxiliary_buses(net, ppc, bus_lookup, "trafo3w", "hv_bus", aux) ``` 331 1 ``` net["_pd2ppc_lookups"]["bus"] = bus_lookup ``` 332 1 ``` net["_pd2ppc_lookups"]["aux"] = aux ``` 333 334 335 1 ```def _fill_auxiliary_buses(net, ppc, bus_lookup, element, bus_column, aux): ``` 336 1 ``` element_bus_idx = bus_lookup[net[element][bus_column].values] ``` 337 1 ``` aux_idx = bus_lookup[aux[element]] ``` 338 1 ``` ppc["bus"][aux_idx, BASE_KV] = ppc["bus"][element_bus_idx, BASE_KV] ``` 339 1 ``` if net._options["mode"] == "opf": ``` 340 1 ``` ppc["bus"][aux_idx, VMIN] = ppc["bus"][element_bus_idx, VMIN] ``` 341 1 ``` ppc["bus"][aux_idx, VMAX] = ppc["bus"][element_bus_idx, VMAX] ``` 342 343 1 ``` if net._options["init_vm_pu"] == "results": ``` 344 1 ``` ppc["bus"][aux_idx, VM] = net["res_%s" % element]["vm_internal_pu"].values ``` 345 ``` else: ``` 346 1 ``` ppc["bus"][aux_idx, VM] = ppc["bus"][element_bus_idx, VM] ``` 347 1 ``` if net._options["init_va_degree"] == "results": ``` 348 1 ``` ppc["bus"][aux_idx, VA] = net["res_%s" % element]["va_internal_degree"].values ``` 349 ``` else: ``` 350 1 ``` ppc["bus"][aux_idx, VA] = ppc["bus"][element_bus_idx, VA] ``` 351 352 353 1 ```def set_reference_buses(net, ppc, bus_lookup, mode): ``` 354 1 ``` if mode == "nx": ``` 355 0 ``` return ``` 356 1 ``` eg_buses = bus_lookup[net.ext_grid.bus.values[net._is_elements["ext_grid"]]] ``` 357 1 ``` ppc["bus"][eg_buses, BUS_TYPE] = REF ``` 358 1 ``` if mode == "sc": ``` 359 1 ``` gen_slacks = net._is_elements["gen"] # generators are slacks for short-circuit calculation ``` 360 ``` else: ``` 361 1 ``` gen_slacks = net._is_elements["gen"] & net.gen["slack"].values ``` 362 1 ``` if gen_slacks.any(): ``` 363 1 ``` slack_buses = net.gen["bus"].values[gen_slacks] ``` 364 1 ``` ppc["bus"][bus_lookup[slack_buses], BUS_TYPE] = REF ``` 365 366 367 1 ```def _calc_pq_elements_and_add_on_ppc(net, ppc, sequence= None): ``` 368 ``` # init values ``` 369 1 ``` b, p, q = np.array([], dtype=int), np.array([]), np.array([]) ``` 370 371 1 ``` _is_elements = net["_is_elements"] ``` 372 1 ``` voltage_depend_loads = net["_options"]["voltage_depend_loads"] ``` 373 1 ``` mode = net["_options"]["mode"] ``` 374 1 ``` pq_elements = ["load", "motor", "sgen", "storage", "ward", "xward"] ``` 375 1 ``` for element in pq_elements: ``` 376 1 ``` tab = net[element] ``` 377 1 ``` if len(tab): ``` 378 1 ``` if element == "load" and voltage_depend_loads: ``` 379 1 ``` cz = tab["const_z_percent"].values / 100. ``` 380 1 ``` ci = tab["const_i_percent"].values / 100. ``` 381 1 ``` if ((cz + ci) > 1).any(): ``` 382 0 ``` raise ValueError("const_z_percent + const_i_percent need to be less or equal to " + ``` 383 ``` "100%!") ``` 384 385 ``` # cumulative sum of constant-current loads ``` 386 1 ``` b_zip = tab["bus"].values ``` 387 1 ``` load_counter = Counter(b_zip) ``` 388 389 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 390 1 ``` b_zip = bus_lookup[b_zip] ``` 391 1 ``` load_counter = {bus_lookup[k]: v for k, v in load_counter.items()} ``` 392 1 ``` b_zip, ci_sum, cz_sum = _sum_by_group(b_zip, ci, cz) ``` 393 394 1 ``` for bus, no_loads in load_counter.items(): ``` 395 1 ``` ci_sum[b_zip == bus] /= no_loads ``` 396 1 ``` cz_sum[b_zip == bus] /= no_loads ``` 397 398 1 ``` ppc["bus"][b_zip, CID] = ci_sum ``` 399 1 ``` ppc["bus"][b_zip, CZD] = cz_sum ``` 400 1 ``` active = _is_elements[element] ``` 401 1 ``` sign = -1 if element == "sgen" else 1 ``` 402 1 ``` if element == "motor": ``` 403 1 ``` p_mw, q_mvar = _get_motor_pq(net) ``` 404 1 ``` p = np.hstack([p, p_mw]) ``` 405 1 ``` q = np.hstack([q, q_mvar]) ``` 406 1 ``` elif element.endswith("ward"): ``` 407 1 ``` p = np.hstack([p, tab["ps_mw"].values * active * sign]) ``` 408 1 ``` q = np.hstack([q, tab["qs_mvar"].values * active * sign]) ``` 409 ``` else: ``` 410 1 ``` scaling = tab["scaling"].values ``` 411 1 ``` p = np.hstack([p, tab["p_mw"].values * active * scaling * sign]) ``` 412 1 ``` q = np.hstack([q, tab["q_mvar"].values * active * scaling * sign]) ``` 413 1 ``` b = np.hstack([b, tab["bus"].values]) ``` 414 415 1 ``` for element in ["asymmetric_load", "asymmetric_sgen"]: ``` 416 1 ``` if len(net[element]) > 0 and mode == "pf": ``` 417 1 ``` p_mw, q_mvar = _get_symmetric_pq_of_unsymetric_element(net, element) ``` 418 1 ``` sign = -1 if element.endswith("sgen") else 1 ``` 419 1 ``` p = np.hstack([p, p_mw * sign]) ``` 420 1 ``` q = np.hstack([q, q_mvar * sign ]) ``` 421 1 ``` b = np.hstack([b, net[element]["bus"].values]) ``` 422 423 ``` # sum up p & q of bus elements ``` 424 1 ``` if b.size: ``` 425 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 426 1 ``` b = bus_lookup[b] ``` 427 1 ``` b, vp, vq = _sum_by_group(b, p, q) ``` 428 1 ``` ppc["bus"][b, PD] = vp ``` 429 1 ``` ppc["bus"][b, QD] = vq ``` 430 431 432 1 ```def _get_symmetric_pq_of_unsymetric_element(net, element): ``` 433 1 ``` scale = net["_is_elements"][element] * net[element]["scaling"].values.T ``` 434 1 ``` q_mvar = np.sum(net[element][["q_a_mvar", "q_b_mvar", "q_c_mvar"]].values, axis=1) ``` 435 1 ``` p_mw = np.sum(net[element][["p_a_mw", "p_b_mw", "p_c_mw"]].values, axis=1) ``` 436 1 ``` return p_mw*scale, q_mvar*scale ``` 437 438 1 ```def _get_motor_pq(net): ``` 439 1 ``` tab = net["motor"] ``` 440 1 ``` active = net._is_elements["motor"] ``` 441 1 ``` scale = tab["loading_percent"].values/100 *tab["scaling"].values*active ``` 442 443 1 ``` efficiency = tab["efficiency_percent"].values ``` 444 1 ``` p_mech = tab["pn_mech_mw"].values ``` 445 1 ``` cos_phi = tab["cos_phi"].values ``` 446 447 1 ``` p_mw = p_mech / efficiency * 100 * scale ``` 448 1 ``` s_mvar = p_mw / cos_phi ``` 449 1 ``` q_mvar = np.sqrt(s_mvar**2 - p_mw**2) ``` 450 1 ``` return p_mw, q_mvar ``` 451 452 1 ```def _calc_shunts_and_add_on_ppc(net, ppc): ``` 453 ``` # init values ``` 454 1 ``` b, p, q = np.array([], dtype=int), np.array([]), np.array([]) ``` 455 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 456 ``` # get in service elements ``` 457 1 ``` _is_elements = net["_is_elements"] ``` 458 459 1 ``` s = net["shunt"] ``` 460 1 ``` if len(s) > 0: ``` 461 1 ``` vl = _is_elements["shunt"] ``` 462 1 ``` v_ratio = (ppc["bus"][bus_lookup[s["bus"].values], BASE_KV] / s["vn_kv"].values) ** 2 ``` 463 1 ``` q = np.hstack([q, s["q_mvar"].values * s["step"].values * v_ratio * vl]) ``` 464 1 ``` p = np.hstack([p, s["p_mw"].values * s["step"].values * v_ratio * vl]) ``` 465 1 ``` b = np.hstack([b, s["bus"].values]) ``` 466 467 1 ``` w = net["ward"] ``` 468 1 ``` if len(w) > 0: ``` 469 1 ``` vl = _is_elements["ward"] ``` 470 1 ``` q = np.hstack([q, w["qz_mvar"].values * vl]) ``` 471 1 ``` p = np.hstack([p, w["pz_mw"].values * vl]) ``` 472 1 ``` b = np.hstack([b, w["bus"].values]) ``` 473 474 1 ``` xw = net["xward"] ``` 475 1 ``` if len(xw) > 0: ``` 476 1 ``` vl = _is_elements["xward"] ``` 477 1 ``` q = np.hstack([q, xw["qz_mvar"].values * vl]) ``` 478 1 ``` p = np.hstack([p, xw["pz_mw"].values * vl]) ``` 479 1 ``` b = np.hstack([b, xw["bus"].values]) ``` 480 481 1 ``` loss_location = net._options["trafo3w_losses"].lower() ``` 482 1 ``` trafo3w = net["trafo3w"] ``` 483 1 ``` if loss_location == "star" and len(trafo3w) > 0: ``` 484 1 ``` pfe_mw = trafo3w["pfe_kw"].values * 1e-3 ``` 485 1 ``` i0 = trafo3w["i0_percent"].values ``` 486 1 ``` sn_mva = trafo3w["sn_hv_mva"].values ``` 487 488 1 ``` q_mvar = (sn_mva * i0 / 100.) ** 2 - pfe_mw ** 2 ``` 489 1 ``` q_mvar[q_mvar < 0] = 0 ``` 490 1 ``` q_mvar = np.sqrt(q_mvar) ``` 491 492 1 ``` vn_hv_trafo = trafo3w["vn_hv_kv"].values ``` 493 1 ``` vn_hv_bus = ppc["bus"][bus_lookup[trafo3w.hv_bus.values], BASE_KV] ``` 494 1 ``` v_ratio = (vn_hv_bus / vn_hv_trafo) ** 2 ``` 495 496 1 ``` q = np.hstack([q, q_mvar * v_ratio]) ``` 497 1 ``` p = np.hstack([p, pfe_mw * v_ratio]) ``` 498 1 ``` b = np.hstack([b, net._pd2ppc_lookups["aux"]["trafo3w"]]) ``` 499 500 ``` # if array is not empty ``` 501 1 ``` if b.size: ``` 502 1 ``` b = bus_lookup[b] ``` 503 1 ``` b, vp, vq = _sum_by_group(b, p, q) ``` 504 505 1 ``` ppc["bus"][b, GS] = vp ``` 506 1 ``` ppc["bus"][b, BS] = -vq ``` 507 508 509 1 ```def _add_gen_impedances_ppc(net, ppc): ``` 510 1 ``` _add_ext_grid_sc_impedance(net, ppc) ``` 511 1 ``` _add_gen_sc_impedance(net, ppc) ``` 512 513 514 1 ```def _add_ext_grid_sc_impedance(net, ppc): ``` 515 1 ``` from pandapower.shortcircuit.idx_bus import C_MAX, C_MIN ``` 516 1 ``` mode = net._options["mode"] ``` 517 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 518 1 ``` if mode == "sc": ``` 519 1 ``` case = net._options["case"] ``` 520 ``` else: ``` 521 1 ``` case = "max" ``` 522 1 ``` eg = net["ext_grid"][net._is_elements["ext_grid"]] ``` 523 1 ``` if len(eg) == 0: ``` 524 1 ``` return ``` 525 1 ``` eg_buses = eg.bus.values ``` 526 1 ``` eg_buses_ppc = bus_lookup[eg_buses] ``` 527 528 1 ``` if mode == "sc": ``` 529 1 ``` c = ppc["bus"][eg_buses_ppc, C_MAX] if case == "max" else ppc["bus"][eg_buses_ppc, C_MIN] ``` 530 ``` else: ``` 531 1 ``` c = 1.1 ``` 532 1 ``` if not "s_sc_%s_mva" % case in eg: ``` 533 0 ``` raise ValueError("short circuit apparent power s_sc_%s_mva needs to be specified for "% case + ``` 534 ``` "external grid \n Try: net.ext_grid['s_sc_max_mva'] = 1000" ) ``` 535 1 ``` s_sc = eg["s_sc_%s_mva" % case].values/ppc['baseMVA'] ``` 536 1 ``` if not "rx_%s" % case in eg: ``` 537 0 ``` raise ValueError("short circuit R/X rate rx_%s needs to be specified for external grid \n Try: net.ext_grid['rx_max'] = 0.1" % ``` 538 ``` case) ``` 539 1 ``` rx = eg["rx_%s" % case].values ``` 540 541 1 ``` z_grid = c / s_sc ``` 542 1 ``` if mode == 'pf_3ph': ``` 543 1 ``` z_grid = c / (s_sc/3) # 3 phase power divided to get 1 ph power ``` 544 1 ``` x_grid = z_grid / np.sqrt(rx ** 2 + 1) ``` 545 1 ``` r_grid = rx * x_grid ``` 546 1 ``` eg["r"] = r_grid ``` 547 1 ``` eg["x"] = x_grid ``` 548 549 1 ``` y_grid = 1 / (r_grid + x_grid * 1j) ``` 550 1 ``` buses, gs, bs = _sum_by_group(eg_buses_ppc, y_grid.real, y_grid.imag) ``` 551 1 ``` ppc["bus"][buses, GS] = gs * ppc['baseMVA'] ``` 552 1 ``` ppc["bus"][buses, BS] = bs * ppc['baseMVA'] ``` 553 1 ``` return gs * ppc['baseMVA'], bs * ppc['baseMVA'] ``` 554 555 556 1 ```def _add_gen_sc_impedance(net, ppc): ``` 557 1 ``` from pandapower.shortcircuit.idx_bus import C_MAX ``` 558 1 ``` gen = net["gen"][net._is_elements["gen"]] ``` 559 1 ``` if len(gen) == 0: ``` 560 1 ``` return ``` 561 1 ``` gen_buses = gen.bus.values ``` 562 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 563 1 ``` gen_buses_ppc = bus_lookup[gen_buses] ``` 564 1 ``` vn_gen = gen.vn_kv.values ``` 565 1 ``` sn_gen = gen.sn_mva.values ``` 566 567 1 ``` rdss_pu = gen.rdss_pu.values ``` 568 1 ``` xdss_pu = gen.xdss_pu.values ``` 569 1 ``` gens_without_r = np.isnan(rdss_pu) ``` 570 1 ``` if gens_without_r.any(): ``` 571 ``` # use the estimations from the IEC standard for generators without defined rdss_pu ``` 572 1 ``` lv_gens = (vn_gen <= 1.) & gens_without_r ``` 573 1 ``` hv_gens = (vn_gen > 1.) & gens_without_r ``` 574 1 ``` large_hv_gens = (sn_gen >= 100) & hv_gens ``` 575 1 ``` small_hv_gens = (sn_gen < 100) & hv_gens ``` 576 1 ``` rdss_pu[lv_gens] = 0.15 * xdss_pu[lv_gens] ``` 577 1 ``` rdss_pu[large_hv_gens] = 0.05 * xdss_pu[large_hv_gens] ``` 578 1 ``` rdss_pu[small_hv_gens] = 0.07 * xdss_pu[small_hv_gens] ``` 579 580 1 ``` vn_net = ppc["bus"][gen_buses_ppc, BASE_KV] ``` 581 1 ``` cmax = ppc["bus"][gen_buses_ppc, C_MAX] ``` 582 1 ``` phi_gen = np.arccos(gen.cos_phi.values) ``` 583 1 ``` kg = vn_gen / vn_net * cmax / (1 + xdss_pu * np.sin(phi_gen)) ``` 584 585 1 ``` z_gen = (rdss_pu + xdss_pu * 1j) * kg / sn_gen ``` 586 1 ``` y_gen = 1 / z_gen ``` 587 588 1 ``` buses, gs, bs = _sum_by_group(gen_buses_ppc, y_gen.real, y_gen.imag) ``` 589 1 ``` ppc["bus"][buses, GS] = gs ``` 590 1 ``` ppc["bus"][buses, BS] = bs ``` 591 592 593 1 ```def _add_motor_impedances_ppc(net, ppc): ``` 594 1 ``` if net._options["case"] == "min": ``` 595 1 ``` return ``` 596 1 ``` motor = net["motor"][net._is_elements["motor"]] ``` 597 1 ``` if motor.empty: ``` 598 1 ``` return ``` 599 1 ``` for par in ["vn_kv", "lrc_pu", "efficiency_n_percent", "cos_phi_n", "rx", "pn_mech_mw"]: ``` 600 1 ``` if any(pd.isnull(motor[par])): ``` 601 0 ``` raise UserWarning("%s needs to be specified for all motors in net.motor.%s" % (par, par)) ``` 602 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 603 1 ``` motor_buses_ppc = bus_lookup[motor.bus.values] ``` 604 1 ``` vn_net = ppc["bus"][motor_buses_ppc, BASE_KV] ``` 605 606 1 ``` efficiency = motor.efficiency_n_percent.values ``` 607 1 ``` cos_phi = motor.cos_phi_n.values ``` 608 1 ``` p_mech = motor.pn_mech_mw.values ``` 609 1 ``` vn_kv = motor.vn_kv.values ``` 610 1 ``` lrc = motor.lrc_pu.values ``` 611 1 ``` rx = motor.rx.values ``` 612 613 1 ``` s_motor = p_mech / (efficiency/100 * cos_phi) ``` 614 1 ``` z_motor_ohm = 1 / lrc * vn_kv**2 / s_motor ``` 615 1 ``` z_motor_pu = z_motor_ohm / (vn_net**2 / net.sn_mva) ``` 616 617 1 ``` x_motor_pu = z_motor_pu / np.sqrt(rx ** 2 + 1) ``` 618 1 ``` r_motor_pu = rx * x_motor_pu ``` 619 1 ``` y_motor_pu = 1 / (r_motor_pu + x_motor_pu * 1j) ``` 620 621 1 ``` buses, gs, bs = _sum_by_group(motor_buses_ppc, y_motor_pu.real, y_motor_pu.imag) ``` 622 1 ``` ppc["bus"][buses, GS] = gs ``` 623 1 ``` ppc["bus"][buses, BS] = bs ``` 624 625 626 1 ```def _add_c_to_ppc(net, ppc): ``` 627 1 ``` from pandapower.shortcircuit.idx_bus import C_MAX, C_MIN ``` 628 1 ``` ppc["bus"][:, C_MAX] = 1.1 ``` 629 1 ``` ppc["bus"][:, C_MIN] = 1. ``` 630 1 ``` lv_buses = np.where(ppc["bus"][:, BASE_KV] < 1.) ``` 631 1 ``` if len(lv_buses) > 0: ``` 632 1 ``` lv_tol_percent = net["_options"]["lv_tol_percent"] ``` 633 1 ``` if lv_tol_percent == 10: ``` 634 1 ``` c_ns = 1.1 ``` 635 1 ``` elif lv_tol_percent == 6: ``` 636 1 ``` c_ns = 1.05 ``` 637 ``` else: ``` 638 0 ``` raise ValueError("Voltage tolerance in the low voltage grid has" + ``` 639 ``` " to be either 6% or 10% according to IEC 60909") ``` 640 1 ``` ppc["bus"][lv_buses, C_MAX] = c_ns ``` 641 1 ``` ppc["bus"][lv_buses, C_MIN] = .95 ```

Read our documentation on viewing source code .