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, phase_to_sequence ``` 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, sequence=None): ``` 205 1 ``` if isinstance(init_v, str): ``` 206 1 ``` if init_v == "results": ``` 207 1 ``` res_table = "res_bus" if sequence is None else "res_bus_3ph" ``` 208 209 ``` # init voltage possible if bus results are available ``` 210 1 ``` if res_table not in net or not net[res_table].index.equals(net.bus.index): ``` 211 ``` # cannot init from results, since sorting of results is different from element table ``` 212 ``` # TO BE REVIEWED! Why there was no raise before this commit? ``` 213 0 ``` raise UserWarning("Init from results not possible. Index of %s do not match with bus. " ``` 214 ``` "You should sort res_bus before calling runpp." % res_table) ``` 215 216 1 ``` if res_table == "res_bus_3ph": ``` 217 1 ``` vm = net.res_bus_3ph[["vm_a_pu", "vm_b_pu", "vm_c_pu"]].values.T ``` 218 1 ``` va = net.res_bus_3ph[["va_a_degree", "va_b_degree", "va_c_degree"]].values.T ``` 219 220 1 ``` voltage_vector = phase_to_sequence(vm * np.exp(1j * np.pi * va / 180.))[sequence, :] ``` 221 222 1 ``` if mode == "magnitude": ``` 223 1 ``` return np.abs(voltage_vector) ``` 224 1 ``` elif mode == "angle": ``` 225 1 ``` return np.angle(voltage_vector) * 180 / np.pi ``` 226 ``` else: ``` 227 0 ``` raise UserWarning(str(mode)+" for initialization not available!") ``` 228 ``` else: ``` 229 1 ``` if mode == "magnitude": ``` 230 1 ``` return net[res_table]["vm_pu"].values.copy() ``` 231 1 ``` elif mode == "angle": ``` 232 1 ``` return net[res_table]["va_degree"].values.copy() ``` 233 ``` else: ``` 234 0 ``` raise UserWarning(str(mode)+" for initialization not available!") ``` 235 1 ``` if init_v == "flat": ``` 236 1 ``` return None ``` 237 1 ``` elif isinstance(init_v, (float, np.ndarray, list)) and sequence is None or sequence == 1: ``` 238 1 ``` return init_v ``` 239 1 ``` elif isinstance(init_v, pd.Series) and sequence is None or sequence == 1: ``` 240 0 ``` if init_v.index.equals(net.bus.index): ``` 241 0 ``` return init_v.loc[net.bus.index] ``` 242 ``` else: ``` 243 0 ``` raise UserWarning("Voltage starting vector indices do not match bus indices") ``` 244 245 246 1 ```def _build_bus_ppc(net, ppc, sequence=None): ``` 247 ``` """ ``` 248 ``` Generates the ppc["bus"] array and the lookup pandapower indices -> ppc indices ``` 249 ``` """ ``` 250 1 ``` init_vm_pu = net["_options"]["init_vm_pu"] ``` 251 1 ``` init_va_degree = net["_options"]["init_va_degree"] ``` 252 1 ``` mode = net["_options"]["mode"] ``` 253 1 ``` numba = net["_options"]["numba"] if "numba" in net["_options"] else False ``` 254 255 ``` # get bus indices ``` 256 1 ``` nr_xward = len(net.xward) ``` 257 1 ``` nr_trafo3w = len(net.trafo3w) ``` 258 1 ``` aux = dict() ``` 259 1 ``` if nr_xward > 0 or nr_trafo3w > 0: ``` 260 1 ``` bus_indices = [net["bus"].index.values, np.array([], dtype=np.int64)] ``` 261 1 ``` max_idx = max(net["bus"].index) + 1 ``` 262 1 ``` if nr_xward > 0: ``` 263 1 ``` aux_xward = np.arange(max_idx, max_idx + nr_xward, dtype=np.int64) ``` 264 1 ``` aux["xward"] = aux_xward ``` 265 1 ``` bus_indices.append(aux_xward) ``` 266 1 ``` if nr_trafo3w: ``` 267 1 ``` aux_trafo3w = np.arange(max_idx + nr_xward, max_idx + nr_xward + nr_trafo3w) ``` 268 1 ``` aux["trafo3w"] = aux_trafo3w ``` 269 1 ``` bus_indices.append(aux_trafo3w) ``` 270 1 ``` bus_index = np.concatenate(bus_indices) ``` 271 ``` else: ``` 272 1 ``` bus_index = net["bus"].index.values ``` 273 ``` # get in service elements ``` 274 275 1 ``` if mode == "nx": ``` 276 1 ``` bus_lookup = create_consecutive_bus_lookup(net, bus_index) ``` 277 ``` else: ``` 278 1 ``` _is_elements = net["_is_elements"] ``` 279 1 ``` eg_is_mask = _is_elements['ext_grid'] ``` 280 1 ``` gen_is_mask = _is_elements['gen'] ``` 281 1 ``` bus_is_idx = _is_elements['bus_is_idx'] ``` 282 1 ``` bus_lookup = create_bus_lookup(net, bus_index, bus_is_idx, ``` 283 ``` gen_is_mask, eg_is_mask, numba=numba) ``` 284 285 1 ``` n_bus_ppc = len(bus_index) ``` 286 ``` # init ppc with empty values ``` 287 1 ``` ppc["bus"] = np.zeros(shape=(n_bus_ppc, bus_cols), dtype=float) ``` 288 1 ``` ppc["bus"][:, :15] = np.array([0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 2, 0, 0., 0.]) # changes of ``` 289 ``` # voltage limits (2 and 0) must be considered in check_opf_data ``` 290 291 1 ``` if sequence is not None and sequence != 1: ``` 292 1 ``` ppc["bus"][:, VM] = 0. ``` 293 294 1 ``` if mode == "sc": ``` 295 1 ``` from pandapower.shortcircuit.idx_bus import bus_cols_sc ``` 296 1 ``` bus_sc = np.empty(shape=(n_bus_ppc, bus_cols_sc), dtype=float) ``` 297 1 ``` bus_sc.fill(np.nan) ``` 298 1 ``` ppc["bus"] = np.hstack((ppc["bus"], bus_sc)) ``` 299 300 ``` # apply consecutive bus numbers ``` 301 1 ``` ppc["bus"][:, BUS_I] = np.arange(n_bus_ppc) ``` 302 303 1 ``` n_bus = len(net.bus.index) ``` 304 ``` # init voltages from net ``` 305 1 ``` ppc["bus"][:n_bus, BASE_KV] = net["bus"]["vn_kv"].values ``` 306 ``` # set buses out of service (BUS_TYPE == 4) ``` 307 1 ``` if nr_xward > 0 or nr_trafo3w > 0: ``` 308 1 ``` in_service = np.concatenate([net["bus"]["in_service"].values, ``` 309 ``` net["xward"]["in_service"].values, ``` 310 ``` net["trafo3w"]["in_service"].values]) ``` 311 ``` else: ``` 312 1 ``` in_service = net["bus"]["in_service"].values ``` 313 1 ``` ppc["bus"][~in_service, BUS_TYPE] = NONE ``` 314 1 ``` if mode != "nx": ``` 315 1 ``` set_reference_buses(net, ppc, bus_lookup, mode) ``` 316 1 ``` vm_pu = get_voltage_init_vector(net, init_vm_pu, "magnitude", sequence=sequence) ``` 317 1 ``` if vm_pu is not None: ``` 318 1 ``` ppc["bus"][:n_bus, VM] = vm_pu ``` 319 320 1 ``` va_degree = get_voltage_init_vector(net, init_va_degree, "angle", sequence=sequence) ``` 321 1 ``` if va_degree is not None: ``` 322 1 ``` ppc["bus"][:n_bus, VA] = va_degree ``` 323 324 1 ``` if mode == "sc": ``` 325 1 ``` _add_c_to_ppc(net, ppc) ``` 326 327 1 ``` if net._options["mode"] == "opf": ``` 328 1 ``` if "max_vm_pu" in net.bus: ``` 329 1 ``` ppc["bus"][:n_bus, VMAX] = net["bus"].max_vm_pu.values ``` 330 ``` else: ``` 331 1 ``` ppc["bus"][:n_bus, VMAX] = 2. # changes of VMAX must be considered in check_opf_data ``` 332 1 ``` if "min_vm_pu" in net.bus: ``` 333 1 ``` ppc["bus"][:n_bus, VMIN] = net["bus"].min_vm_pu.values ``` 334 ``` else: ``` 335 1 ``` ppc["bus"][:n_bus, VMIN] = 0. # changes of VMIN must be considered in check_opf_data ``` 336 337 1 ``` if len(net.xward): ``` 338 1 ``` _fill_auxiliary_buses(net, ppc, bus_lookup, "xward", "bus", aux) ``` 339 340 1 ``` if len(net.trafo3w): ``` 341 1 ``` _fill_auxiliary_buses(net, ppc, bus_lookup, "trafo3w", "hv_bus", aux) ``` 342 1 ``` net["_pd2ppc_lookups"]["bus"] = bus_lookup ``` 343 1 ``` net["_pd2ppc_lookups"]["aux"] = aux ``` 344 345 346 1 ```def _fill_auxiliary_buses(net, ppc, bus_lookup, element, bus_column, aux): ``` 347 1 ``` element_bus_idx = bus_lookup[net[element][bus_column].values] ``` 348 1 ``` aux_idx = bus_lookup[aux[element]] ``` 349 1 ``` ppc["bus"][aux_idx, BASE_KV] = ppc["bus"][element_bus_idx, BASE_KV] ``` 350 1 ``` if net._options["mode"] == "opf": ``` 351 1 ``` ppc["bus"][aux_idx, VMIN] = ppc["bus"][element_bus_idx, VMIN] ``` 352 1 ``` ppc["bus"][aux_idx, VMAX] = ppc["bus"][element_bus_idx, VMAX] ``` 353 354 1 ``` if net._options["init_vm_pu"] == "results": ``` 355 1 ``` ppc["bus"][aux_idx, VM] = net["res_%s" % element]["vm_internal_pu"].values ``` 356 ``` else: ``` 357 1 ``` ppc["bus"][aux_idx, VM] = ppc["bus"][element_bus_idx, VM] ``` 358 1 ``` if net._options["init_va_degree"] == "results": ``` 359 1 ``` ppc["bus"][aux_idx, VA] = net["res_%s" % element]["va_internal_degree"].values ``` 360 ``` else: ``` 361 1 ``` ppc["bus"][aux_idx, VA] = ppc["bus"][element_bus_idx, VA] ``` 362 363 364 1 ```def set_reference_buses(net, ppc, bus_lookup, mode): ``` 365 1 ``` if mode == "nx": ``` 366 0 ``` return ``` 367 1 ``` eg_buses = bus_lookup[net.ext_grid.bus.values[net._is_elements["ext_grid"]]] ``` 368 1 ``` ppc["bus"][eg_buses, BUS_TYPE] = REF ``` 369 1 ``` if mode == "sc": ``` 370 1 ``` gen_slacks = net._is_elements["gen"] # generators are slacks for short-circuit calculation ``` 371 ``` else: ``` 372 1 ``` gen_slacks = net._is_elements["gen"] & net.gen["slack"].values ``` 373 1 ``` if gen_slacks.any(): ``` 374 1 ``` slack_buses = net.gen["bus"].values[gen_slacks] ``` 375 1 ``` ppc["bus"][bus_lookup[slack_buses], BUS_TYPE] = REF ``` 376 377 378 1 ```def _calc_pq_elements_and_add_on_ppc(net, ppc, sequence= None): ``` 379 ``` # init values ``` 380 1 ``` b, p, q = np.array([], dtype=int), np.array([]), np.array([]) ``` 381 382 1 ``` _is_elements = net["_is_elements"] ``` 383 1 ``` voltage_depend_loads = net["_options"]["voltage_depend_loads"] ``` 384 1 ``` mode = net["_options"]["mode"] ``` 385 1 ``` pq_elements = ["load", "motor", "sgen", "storage", "ward", "xward"] ``` 386 1 ``` for element in pq_elements: ``` 387 1 ``` tab = net[element] ``` 388 1 ``` if len(tab): ``` 389 1 ``` if element == "load" and voltage_depend_loads: ``` 390 1 ``` cz = tab["const_z_percent"].values / 100. ``` 391 1 ``` ci = tab["const_i_percent"].values / 100. ``` 392 1 ``` if ((cz + ci) > 1).any(): ``` 393 0 ``` raise ValueError("const_z_percent + const_i_percent need to be less or equal to " + ``` 394 ``` "100%!") ``` 395 396 ``` # cumulative sum of constant-current loads ``` 397 1 ``` b_zip = tab["bus"].values ``` 398 1 ``` load_counter = Counter(b_zip) ``` 399 400 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 401 1 ``` b_zip = bus_lookup[b_zip] ``` 402 1 ``` load_counter = {bus_lookup[k]: v for k, v in load_counter.items()} ``` 403 1 ``` b_zip, ci_sum, cz_sum = _sum_by_group(b_zip, ci, cz) ``` 404 405 1 ``` for bus, no_loads in load_counter.items(): ``` 406 1 ``` ci_sum[b_zip == bus] /= no_loads ``` 407 1 ``` cz_sum[b_zip == bus] /= no_loads ``` 408 409 1 ``` ppc["bus"][b_zip, CID] = ci_sum ``` 410 1 ``` ppc["bus"][b_zip, CZD] = cz_sum ``` 411 1 ``` active = _is_elements[element] ``` 412 1 ``` sign = -1 if element == "sgen" else 1 ``` 413 1 ``` if element == "motor": ``` 414 1 ``` p_mw, q_mvar = _get_motor_pq(net) ``` 415 1 ``` p = np.hstack([p, p_mw]) ``` 416 1 ``` q = np.hstack([q, q_mvar]) ``` 417 1 ``` elif element.endswith("ward"): ``` 418 1 ``` p = np.hstack([p, tab["ps_mw"].values * active * sign]) ``` 419 1 ``` q = np.hstack([q, tab["qs_mvar"].values * active * sign]) ``` 420 ``` else: ``` 421 1 ``` scaling = tab["scaling"].values ``` 422 1 ``` p = np.hstack([p, tab["p_mw"].values * active * scaling * sign]) ``` 423 1 ``` q = np.hstack([q, tab["q_mvar"].values * active * scaling * sign]) ``` 424 1 ``` b = np.hstack([b, tab["bus"].values]) ``` 425 426 1 ``` for element in ["asymmetric_load", "asymmetric_sgen"]: ``` 427 1 ``` if len(net[element]) > 0 and mode == "pf": ``` 428 1 ``` p_mw, q_mvar = _get_symmetric_pq_of_unsymetric_element(net, element) ``` 429 1 ``` sign = -1 if element.endswith("sgen") else 1 ``` 430 1 ``` p = np.hstack([p, p_mw * sign]) ``` 431 1 ``` q = np.hstack([q, q_mvar * sign ]) ``` 432 1 ``` b = np.hstack([b, net[element]["bus"].values]) ``` 433 434 ``` # sum up p & q of bus elements ``` 435 1 ``` if b.size: ``` 436 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 437 1 ``` b = bus_lookup[b] ``` 438 1 ``` b, vp, vq = _sum_by_group(b, p, q) ``` 439 1 ``` ppc["bus"][b, PD] = vp ``` 440 1 ``` ppc["bus"][b, QD] = vq ``` 441 442 443 1 ```def _get_symmetric_pq_of_unsymetric_element(net, element): ``` 444 1 ``` scale = net["_is_elements"][element] * net[element]["scaling"].values.T ``` 445 1 ``` q_mvar = np.sum(net[element][["q_a_mvar", "q_b_mvar", "q_c_mvar"]].values, axis=1) ``` 446 1 ``` p_mw = np.sum(net[element][["p_a_mw", "p_b_mw", "p_c_mw"]].values, axis=1) ``` 447 1 ``` return p_mw*scale, q_mvar*scale ``` 448 449 1 ```def _get_motor_pq(net): ``` 450 1 ``` tab = net["motor"] ``` 451 1 ``` active = net._is_elements["motor"] ``` 452 1 ``` scale = tab["loading_percent"].values/100 *tab["scaling"].values*active ``` 453 454 1 ``` efficiency = tab["efficiency_percent"].values ``` 455 1 ``` p_mech = tab["pn_mech_mw"].values ``` 456 1 ``` cos_phi = tab["cos_phi"].values ``` 457 458 1 ``` p_mw = p_mech / efficiency * 100 * scale ``` 459 1 ``` s_mvar = p_mw / cos_phi ``` 460 1 ``` q_mvar = np.sqrt(s_mvar**2 - p_mw**2) ``` 461 1 ``` return p_mw, q_mvar ``` 462 463 1 ```def _calc_shunts_and_add_on_ppc(net, ppc): ``` 464 ``` # init values ``` 465 1 ``` b, p, q = np.array([], dtype=int), np.array([]), np.array([]) ``` 466 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 467 ``` # get in service elements ``` 468 1 ``` _is_elements = net["_is_elements"] ``` 469 470 1 ``` s = net["shunt"] ``` 471 1 ``` if len(s) > 0: ``` 472 1 ``` vl = _is_elements["shunt"] ``` 473 1 ``` v_ratio = (ppc["bus"][bus_lookup[s["bus"].values], BASE_KV] / s["vn_kv"].values) ** 2 ``` 474 1 ``` q = np.hstack([q, s["q_mvar"].values * s["step"].values * v_ratio * vl]) ``` 475 1 ``` p = np.hstack([p, s["p_mw"].values * s["step"].values * v_ratio * vl]) ``` 476 1 ``` b = np.hstack([b, s["bus"].values]) ``` 477 478 1 ``` w = net["ward"] ``` 479 1 ``` if len(w) > 0: ``` 480 1 ``` vl = _is_elements["ward"] ``` 481 1 ``` q = np.hstack([q, w["qz_mvar"].values * vl]) ``` 482 1 ``` p = np.hstack([p, w["pz_mw"].values * vl]) ``` 483 1 ``` b = np.hstack([b, w["bus"].values]) ``` 484 485 1 ``` xw = net["xward"] ``` 486 1 ``` if len(xw) > 0: ``` 487 1 ``` vl = _is_elements["xward"] ``` 488 1 ``` q = np.hstack([q, xw["qz_mvar"].values * vl]) ``` 489 1 ``` p = np.hstack([p, xw["pz_mw"].values * vl]) ``` 490 1 ``` b = np.hstack([b, xw["bus"].values]) ``` 491 492 1 ``` loss_location = net._options["trafo3w_losses"].lower() ``` 493 1 ``` trafo3w = net["trafo3w"] ``` 494 1 ``` if loss_location == "star" and len(trafo3w) > 0: ``` 495 1 ``` pfe_mw = trafo3w["pfe_kw"].values * 1e-3 ``` 496 1 ``` i0 = trafo3w["i0_percent"].values ``` 497 1 ``` sn_mva = trafo3w["sn_hv_mva"].values ``` 498 499 1 ``` q_mvar = (sn_mva * i0 / 100.) ** 2 - pfe_mw ** 2 ``` 500 1 ``` q_mvar[q_mvar < 0] = 0 ``` 501 1 ``` q_mvar = np.sqrt(q_mvar) ``` 502 503 1 ``` vn_hv_trafo = trafo3w["vn_hv_kv"].values ``` 504 1 ``` vn_hv_bus = ppc["bus"][bus_lookup[trafo3w.hv_bus.values], BASE_KV] ``` 505 1 ``` v_ratio = (vn_hv_bus / vn_hv_trafo) ** 2 ``` 506 507 1 ``` q = np.hstack([q, q_mvar * v_ratio]) ``` 508 1 ``` p = np.hstack([p, pfe_mw * v_ratio]) ``` 509 1 ``` b = np.hstack([b, net._pd2ppc_lookups["aux"]["trafo3w"]]) ``` 510 511 ``` # if array is not empty ``` 512 1 ``` if b.size: ``` 513 1 ``` b = bus_lookup[b] ``` 514 1 ``` b, vp, vq = _sum_by_group(b, p, q) ``` 515 516 1 ``` ppc["bus"][b, GS] = vp ``` 517 1 ``` ppc["bus"][b, BS] = -vq ``` 518 519 ```# Short circuit relevant routines ``` 520 1 ```def _add_ext_grid_sc_impedance(net, ppc): ``` 521 1 ``` from pandapower.shortcircuit.idx_bus import C_MAX, C_MIN ``` 522 1 ``` mode = net._options["mode"] ``` 523 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 524 1 ``` if mode == "sc": ``` 525 1 ``` case = net._options["case"] ``` 526 ``` else: ``` 527 1 ``` case = "max" ``` 528 1 ``` eg = net["ext_grid"][net._is_elements["ext_grid"]] ``` 529 1 ``` if len(eg) == 0: ``` 530 1 ``` return ``` 531 1 ``` eg_buses = eg.bus.values ``` 532 1 ``` eg_buses_ppc = bus_lookup[eg_buses] ``` 533 534 1 ``` if mode == "sc": ``` 535 1 ``` c = ppc["bus"][eg_buses_ppc, C_MAX] if case == "max" else ppc["bus"][eg_buses_ppc, C_MIN] ``` 536 ``` else: ``` 537 1 ``` c = 1.1 ``` 538 1 ``` if not "s_sc_%s_mva" % case in eg: ``` 539 0 ``` raise ValueError("short circuit apparent power s_sc_%s_mva needs to be specified for "% case + ``` 540 ``` "external grid \n Try: net.ext_grid['s_sc_max_mva'] = 1000" ) ``` 541 1 ``` s_sc = eg["s_sc_%s_mva" % case].values/ppc['baseMVA'] ``` 542 1 ``` if not "rx_%s" % case in eg: ``` 543 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" % ``` 544 ``` case) ``` 545 1 ``` rx = eg["rx_%s" % case].values ``` 546 547 1 ``` z_grid = c / s_sc ``` 548 1 ``` if mode == 'pf_3ph': ``` 549 1 ``` z_grid = c / (s_sc/3) # 3 phase power divided to get 1 ph power ``` 550 1 ``` x_grid = z_grid / np.sqrt(rx ** 2 + 1) ``` 551 1 ``` r_grid = rx * x_grid ``` 552 1 ``` eg["r"] = r_grid ``` 553 1 ``` eg["x"] = x_grid ``` 554 555 1 ``` y_grid = 1 / (r_grid + x_grid * 1j) ``` 556 1 ``` buses, gs, bs = _sum_by_group(eg_buses_ppc, y_grid.real, y_grid.imag) ``` 557 1 ``` ppc["bus"][buses, GS] = gs * ppc['baseMVA'] ``` 558 1 ``` ppc["bus"][buses, BS] = bs * ppc['baseMVA'] ``` 559 1 ``` return gs * ppc['baseMVA'], bs * ppc['baseMVA'] ``` 560 561 1 ```def _add_motor_impedances_ppc(net, ppc): ``` 562 1 ``` if net._options["case"] == "min": ``` 563 1 ``` return ``` 564 1 ``` motor = net["motor"][net._is_elements["motor"]] ``` 565 1 ``` if motor.empty: ``` 566 1 ``` return ``` 567 1 ``` for par in ["vn_kv", "lrc_pu", "efficiency_n_percent", "cos_phi_n", "rx", "pn_mech_mw"]: ``` 568 1 ``` if any(pd.isnull(motor[par])): ``` 569 0 ``` raise UserWarning("%s needs to be specified for all motors in net.motor.%s" % (par, par)) ``` 570 1 ``` bus_lookup = net["_pd2ppc_lookups"]["bus"] ``` 571 1 ``` motor_buses_ppc = bus_lookup[motor.bus.values] ``` 572 1 ``` vn_net = ppc["bus"][motor_buses_ppc, BASE_KV] ``` 573 574 1 ``` efficiency = motor.efficiency_n_percent.values ``` 575 1 ``` cos_phi = motor.cos_phi_n.values ``` 576 1 ``` p_mech = motor.pn_mech_mw.values ``` 577 1 ``` vn_kv = motor.vn_kv.values ``` 578 1 ``` lrc = motor.lrc_pu.values ``` 579 1 ``` rx = motor.rx.values ``` 580 581 1 ``` s_motor = p_mech / (efficiency/100 * cos_phi) ``` 582 1 ``` z_motor_ohm = 1 / lrc * vn_kv**2 / s_motor ``` 583 1 ``` z_motor_pu = z_motor_ohm / vn_net**2 ``` 584 585 1 ``` x_motor_pu = z_motor_pu / np.sqrt(rx ** 2 + 1) ``` 586 1 ``` r_motor_pu = rx * x_motor_pu ``` 587 1 ``` y_motor_pu = 1 / (r_motor_pu + x_motor_pu * 1j) ``` 588 589 1 ``` buses, gs, bs = _sum_by_group(motor_buses_ppc, y_motor_pu.real, y_motor_pu.imag) ``` 590 1 ``` ppc["bus"][buses, GS] = gs ``` 591 1 ``` ppc["bus"][buses, BS] = bs ``` 592 593 594 1 ```def _add_c_to_ppc(net, ppc): ``` 595 1 ``` from pandapower.shortcircuit.idx_bus import C_MAX, C_MIN ``` 596 1 ``` ppc["bus"][:, C_MAX] = 1.1 ``` 597 1 ``` ppc["bus"][:, C_MIN] = 1. ``` 598 1 ``` lv_buses = np.where(ppc["bus"][:, BASE_KV] < 1.) ``` 599 1 ``` if len(lv_buses) > 0: ``` 600 1 ``` lv_tol_percent = net["_options"]["lv_tol_percent"] ``` 601 1 ``` if lv_tol_percent == 10: ``` 602 1 ``` c_ns = 1.1 ``` 603 1 ``` elif lv_tol_percent == 6: ``` 604 1 ``` c_ns = 1.05 ``` 605 ``` else: ``` 606 0 ``` raise ValueError("Voltage tolerance in the low voltage grid has" + ``` 607 ``` " to be either 6% or 10% according to IEC 60909") ``` 608 1 ``` ppc["bus"][lv_buses, C_MAX] = c_ns ``` 609 1 ``` ppc["bus"][lv_buses, C_MIN] = .95 ```

Read our documentation on viewing source code .