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
# Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6

7 1
import copy
8 1
import math
9 1
from functools import partial
10

11 1
import numpy as np
12 1
import pandas as pd
13

14 1
from pandapower.auxiliary import get_values
15 1
from pandapower.pypower.idx_brch import F_BUS, T_BUS, BR_R, BR_X, BR_B, TAP, SHIFT, BR_STATUS, RATE_A, \
16
    BR_R_ASYM, BR_X_ASYM, branch_cols
17 1
from pandapower.pypower.idx_bus import BASE_KV, VM, VA
18

19

20 1
def _build_branch_ppc(net, ppc):
21
    """
22
    Takes the empty ppc network and fills it with the branch values. The branch
23
    datatype will be np.complex 128 afterwards.
24

25
    .. note:: The order of branches in the ppc is:
26
            1. Lines
27
            2. Transformers
28
            3. 3W Transformers (each 3W Transformer takes up three branches)
29
            4. Impedances
30
            5. Internal branch for extended ward
31

32
    **INPUT**:
33
        **net** -The pandapower format network
34

35
        **ppc** - The PYPOWER format network to fill in values
36

37
    """
38 1
    length = _initialize_branch_lookup(net)
39 1
    lookup = net._pd2ppc_lookups["branch"]
40 1
    mode = net._options["mode"]
41 1
    ppc["branch"] = np.zeros(shape=(length, branch_cols), dtype=np.complex128)
42 1
    if mode == "sc":
43 1
        from pandapower.shortcircuit.idx_brch import branch_cols_sc
44 1
        branch_sc = np.empty(shape=(length, branch_cols_sc), dtype=float)
45 1
        branch_sc.fill(np.nan)
46 1
        ppc["branch"] = np.hstack((ppc["branch"], branch_sc))
47 1
    ppc["branch"][:, :13] = np.array([0, 0, 0, 0, 0, 250, 250, 250, 1, 0, 1, -360, 360])
48 1
    if "line" in lookup:
49 1
        _calc_line_parameter(net, ppc)
50 1
    if "trafo" in lookup:
51 1
        _calc_trafo_parameter(net, ppc)
52 1
    if "trafo3w" in lookup:
53 1
        _calc_trafo3w_parameter(net, ppc)
54 1
    if "impedance" in lookup:
55 1
        _calc_impedance_parameter(net, ppc)
56 1
    if "xward" in lookup:
57 1
        _calc_xward_parameter(net, ppc)
58 1
    if "switch" in lookup:
59 1
        _calc_switch_parameter(net, ppc)
60

61

62 1
def _initialize_branch_lookup(net):
63 1
    start = 0
64 1
    end = 0
65 1
    net._pd2ppc_lookups["branch"] = {}
66 1
    for element in ["line", "trafo", "trafo3w", "impedance", "xward"]:
67 1
        if len(net[element]) > 0:
68 1
            if element == "trafo3w":
69 1
                end = start + len(net[element]) * 3
70
            else:
71 1
                end = start + len(net[element])
72 1
            net._pd2ppc_lookups["branch"][element] = (start, end)
73 1
            start = end
74 1
    if "_impedance_bb_switches" in net and net._impedance_bb_switches.any():
75 1
        end = start + net._impedance_bb_switches.sum()
76 1
        net._pd2ppc_lookups["branch"]["switch"] = (start, end)
77 1
    return end
78

79

80 1
def _calc_trafo3w_parameter(net, ppc):
81 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
82 1
    branch = ppc["branch"]
83 1
    f, t = net["_pd2ppc_lookups"]["branch"]["trafo3w"]
84 1
    trafo_df = _trafo_df_from_trafo3w(net)
85 1
    hv_bus = get_trafo_values(trafo_df, "hv_bus").astype(int)
86 1
    lv_bus = get_trafo_values(trafo_df, "lv_bus").astype(int)
87 1
    in_service = get_trafo_values(trafo_df, "in_service").astype(int)
88 1
    branch[f:t, F_BUS] = bus_lookup[hv_bus]
89 1
    branch[f:t, T_BUS] = bus_lookup[lv_bus]
90 1
    r, x, y, ratio, shift = _calc_branch_values_from_trafo_df(net, ppc, trafo_df)
91 1
    branch[f:t, BR_R] = r
92 1
    branch[f:t, BR_X] = x
93 1
    branch[f:t, BR_B] = y
94 1
    branch[f:t, TAP] = ratio
95 1
    branch[f:t, SHIFT] = shift
96 1
    branch[f:t, BR_STATUS] = in_service
97 1
    if net["_options"]["mode"] == "opf":
98 1
        if "max_loading_percent" in trafo_df:
99 1
            max_load = get_trafo_values(trafo_df, "max_loading_percent")
100 1
            sn_mva = get_trafo_values(trafo_df, "sn_mva")
101 1
            branch[f:t, RATE_A] = max_load / 100. * sn_mva
102
        else:
103 1
            branch[f:t, RATE_A] = np.nan
104

105

106 1
def _calc_line_parameter(net, ppc, elm="line", ppc_elm="branch"):
107
    """
108
    calculates the line parameter in per unit.
109

110
    **INPUT**:
111
        **net** - The pandapower format network
112

113
        **ppc** - the ppc array
114

115
    **OPTIONAL**:
116
        **elm** - The pandapower element (normally "line")
117

118
        **ppc_elm** - The ppc element (normally "branch")
119

120
    **RETURN**:
121
        **t** - Temporary line parameter. Which is a complex128
122
                Nunmpy array. with the following order:
123
                0:bus_a; 1:bus_b; 2:r_pu; 3:x_pu; 4:b_pu
124
    """
125 1
    f, t = net._pd2ppc_lookups[ppc_elm][elm]
126 1
    branch = ppc[ppc_elm]
127 1
    mode = net["_options"]["mode"]
128 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
129 1
    line = net[elm]
130 1
    from_bus = bus_lookup[line["from_bus"].values]
131 1
    to_bus = bus_lookup[line["to_bus"].values]
132 1
    length_km = line["length_km"].values
133 1
    parallel = line["parallel"].values
134 1
    base_kv = ppc["bus"][from_bus, BASE_KV]
135 1
    baseR = np.square(base_kv) / (3 * net.sn_mva) if mode == "pf_3ph" else   np.square(base_kv) / net.sn_mva 
136 1
    branch[f:t, F_BUS] = from_bus
137 1
    branch[f:t, T_BUS] = to_bus
138 1
    branch[f:t, BR_R] = line["r_ohm_per_km"].values * length_km / baseR / parallel
139 1
    branch[f:t, BR_X] = line["x_ohm_per_km"].values * length_km / baseR / parallel
140

141 1
    if mode == "sc":
142
        # temperature correction
143 1
        if net["_options"]["case"] == "min":
144 1
            branch[f:t, BR_R] *= _end_temperature_correction_factor(net, short_circuit=True)
145
    else:
146
        # temperature correction
147 1
        if net["_options"]["consider_line_temperature"]:
148 1
            branch[f:t, BR_R] *= _end_temperature_correction_factor(net)
149

150 1
        b = 2 * net.f_hz * math.pi * line["c_nf_per_km"].values * 1e-9 * baseR * length_km * parallel
151 1
        g = line["g_us_per_km"].values * 1e-6 * baseR * length_km * parallel
152 1
        branch[f:t, BR_B] = b - g * 1j
153
    # in service of lines
154 1
    branch[f:t, BR_STATUS] = line["in_service"].values
155 1
    if net._options["mode"] == "opf":
156
        # RATE_A is conisdered by the (PowerModels) OPF. If zero -> unlimited
157 1
        max_load = line.max_loading_percent.values if "max_loading_percent" in line else 0.
158 1
        vr = net.bus.loc[line["from_bus"].values, "vn_kv"].values * np.sqrt(3.)
159 1
        max_i_ka = line.max_i_ka.values
160 1
        df = line.df.values
161
        # This calculates the maximum apparent power at 1.0 p.u.
162 1
        branch[f:t, RATE_A] = max_load / 100. * max_i_ka * df * parallel * vr
163

164

165 1
def _calc_trafo_parameter(net, ppc):
166
    '''
167
    Calculates the transformer parameter in per unit.
168

169
    **INPUT**:
170
        **net** - The pandapower format network
171

172
    **RETURN**:
173
        **temp_para** -
174
        Temporary transformer parameter. Which is a np.complex128
175
        Numpy array. with the following order:
176
        0:hv_bus; 1:lv_bus; 2:r_pu; 3:x_pu; 4:b_pu; 5:tab, 6:shift
177
    '''
178

179 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
180 1
    f, t = net["_pd2ppc_lookups"]["branch"]["trafo"]
181 1
    branch = ppc["branch"]
182 1
    trafo = net["trafo"]
183 1
    parallel = trafo["parallel"].values
184 1
    branch[f:t, F_BUS] = bus_lookup[trafo["hv_bus"].values]
185 1
    branch[f:t, T_BUS] = bus_lookup[trafo["lv_bus"].values]
186 1
    r, x, y, ratio, shift = _calc_branch_values_from_trafo_df(net, ppc)
187 1
    branch[f:t, BR_R] = r
188 1
    branch[f:t, BR_X] = x
189 1
    branch[f:t, BR_B] = y
190 1
    branch[f:t, TAP] = ratio
191 1
    branch[f:t, SHIFT] = shift
192 1
    branch[f:t, BR_STATUS] = trafo["in_service"].values
193 1
    if any(trafo.df.values <= 0):
194 1
        raise UserWarning("Rating factor df must be positive. Transformers with false "
195
                          "rating factors: %s" % trafo.query('df<=0').index.tolist())
196 1
    if net._options["mode"] == "opf":
197 1
        max_load = trafo.max_loading_percent.values if "max_loading_percent" in trafo else 0
198 1
        sn_mva = trafo.sn_mva.values
199 1
        df = trafo.df.values
200 1
        branch[f:t, RATE_A] = max_load / 100. * sn_mva * df * parallel
201

202

203 1
def get_trafo_values(trafo_df, par):
204 1
    if isinstance(trafo_df, dict):
205 1
        return trafo_df[par]
206
    else:
207 1
        return trafo_df[par].values
208

209

210 1
def _calc_branch_values_from_trafo_df(net, ppc, trafo_df=None):
211
    """
212
    Calculates the MAT/PYPOWER-branch-attributes from the pandapower trafo dataframe.
213

214
    PYPOWER and MATPOWER uses the PI-model to model transformers.
215
    This function calculates the resistance r, reactance x, complex susceptance c and the tap ratio
216
    according to the given parameters.
217

218
    .. warning:: This function returns the subsceptance b as a complex number
219
        **(-img + -re*i)**. MAT/PYPOWER is only intended to calculate the
220
        imaginary part of the subceptance. However, internally c is
221
        multiplied by i. By using subsceptance in this way, it is possible
222
        to consider the ferromagnetic loss of the coil. Which would
223
        otherwise be neglected.
224

225

226
    .. warning:: Tab switches effect calculation as following:
227
        On **high-voltage** side(=1) -> only **tab** gets adapted.
228
        On **low-voltage** side(=2) -> **tab, x, r** get adapted.
229
        This is consistent with Sincal.
230
        The Sincal method in this case is questionable.
231

232

233
    **INPUT**:
234
        **pd_trafo** - The pandapower format Transformer Dataframe.
235
                        The Transformer modell will only readfrom pd_net
236

237
    **RETURN**:
238
        **temp_para** - Temporary transformer parameter. Which is a complex128
239
                        Nunmpy array. with the following order:
240
                        0:r_pu; 1:x_pu; 2:b_pu; 3:tab;
241

242
    """
243 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
244 1
    if trafo_df is None:
245 1
        trafo_df = net["trafo"]
246 1
    lv_bus = get_trafo_values(trafo_df, "lv_bus")
247 1
    vn_lv = ppc["bus"][bus_lookup[lv_bus], BASE_KV]
248
    ### Construct np.array to parse results in ###
249
    # 0:r_pu; 1:x_pu; 2:b_pu; 3:tab;
250 1
    vn_trafo_hv, vn_trafo_lv, shift = _calc_tap_from_dataframe(net, trafo_df)
251 1
    ratio = _calc_nominal_ratio_from_dataframe(ppc, trafo_df, vn_trafo_hv, vn_trafo_lv,
252
                                               bus_lookup)
253 1
    r, x, y = _calc_r_x_y_from_dataframe(net, trafo_df, vn_trafo_lv, vn_lv, ppc)
254 1
    return r, x, y, ratio, shift
255

256

257 1
def _calc_r_x_y_from_dataframe(net, trafo_df, vn_trafo_lv, vn_lv, ppc):
258 1
    mode = net["_options"]["mode"]
259 1
    trafo_model = net["_options"]["trafo_model"]
260

261 1
    r, x = _calc_r_x_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, net.sn_mva)
262 1
    if mode == "sc":
263 1
        y = 0
264 1
        if isinstance(trafo_df, pd.DataFrame):  # 2w trafo is dataframe, 3w trafo is dict
265 1
            from pandapower.shortcircuit.idx_bus import C_MAX
266 1
            bus_lookup = net._pd2ppc_lookups["bus"]
267 1
            cmax = ppc["bus"][bus_lookup[net.trafo.lv_bus.values], C_MAX]
268 1
            kt = _transformer_correction_factor(trafo_df.vk_percent, trafo_df.vkr_percent,
269
                                                trafo_df.sn_mva, cmax)
270 1
            r *= kt
271 1
            x *= kt
272
    else:
273 1
        y = _calc_y_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, net.sn_mva)
274 1
    if trafo_model == "pi":
275 1
        return r, x, y
276 1
    elif trafo_model == "t":
277 1
        return _wye_delta(r, x, y)
278
    else:
279 0
        raise ValueError("Unkonwn Transformer Model %s - valid values ar 'pi' or 't'" % trafo_model)
280

281

282 1
def _wye_delta(r, x, y):
283
    """
284
    20.05.2016 added by Lothar Löwer
285

286
    Calculate transformer Pi-Data based on T-Data
287

288
    """
289 1
    tidx = np.where(y != 0)
290 1
    za_star = (r[tidx] + x[tidx] * 1j) / 2
291 1
    zc_star = -1j / y[tidx]
292 1
    zSum_triangle = za_star * za_star + 2 * za_star * zc_star
293 1
    zab_triangle = zSum_triangle / zc_star
294 1
    zbc_triangle = zSum_triangle / za_star
295 1
    r[tidx] = zab_triangle.real
296 1
    x[tidx] = zab_triangle.imag
297 1
    y[tidx] = -2j / zbc_triangle
298 1
    return r, x, y
299

300

301 1
def _calc_y_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, sn_mva):
302
    """
303
    Calculate the subsceptance y from the transformer dataframe.
304

305
    INPUT:
306

307
        **trafo** (Dataframe) - The dataframe in net.trafo
308
        which contains transformer calculation values.
309

310
    OUTPUT:
311
        **subsceptance** (1d array, np.complex128) - The subsceptance in pu in
312
        the form (-b_img, -b_real)
313
    """
314
   
315 1
    baseR = np.square(vn_lv) / (3*sn_mva) if mode == 'pf_3ph' else np.square(vn_lv) / sn_mva
316 1
    vn_lv_kv = get_trafo_values(trafo_df, "vn_lv_kv")
317 1
    pfe = get_trafo_values(trafo_df, "pfe_kw") * 1e-3
318 1
    parallel = get_trafo_values(trafo_df, "parallel")
319

320
    ### Calculate subsceptance ###
321
   
322 1
    vnl_squared = (vn_lv_kv ** 2)/3 if mode == 'pf_3ph' else  vn_lv_kv **2
323 1
    b_real = pfe / vnl_squared * baseR
324 1
    i0 = get_trafo_values(trafo_df, "i0_percent")
325 1
    sn = get_trafo_values(trafo_df, "sn_mva")
326 1
    b_img = (i0 / 100. * sn) ** 2 - pfe ** 2
327

328 1
    b_img[b_img < 0] = 0
329 1
    b_img = np.sqrt(b_img) * baseR / vnl_squared
330 1
    y = - b_real * 1j - b_img * np.sign(i0)
331 1
    return y / np.square(vn_trafo_lv / vn_lv_kv) * parallel
332

333

334 1
def _calc_tap_from_dataframe(net, trafo_df):
335
    """
336
    Adjust the nominal voltage vnh and vnl to the active tab position "tap_pos".
337
    If "side" is 1 (high-voltage side) the high voltage vnh is adjusted.
338
    If "side" is 2 (low-voltage side) the low voltage vnl is adjusted
339

340
    INPUT:
341
        **net** - The pandapower format network
342

343
        **trafo** (Dataframe) - The dataframe in pd_net["structure"]["trafo"]
344
        which contains transformer calculation values.
345

346
    OUTPUT:
347
        **vn_hv_kv** (1d array, float) - The adusted high voltages
348

349
        **vn_lv_kv** (1d array, float) - The adjusted low voltages
350

351
        **trafo_shift** (1d array, float) - phase shift angle
352

353
    """
354 1
    calculate_voltage_angles = net["_options"]["calculate_voltage_angles"]
355 1
    mode = net["_options"]["mode"]
356 1
    vnh = copy.copy(get_trafo_values(trafo_df, "vn_hv_kv").astype(float))
357 1
    vnl = copy.copy(get_trafo_values(trafo_df, "vn_lv_kv").astype(float))
358 1
    trafo_shift = get_trafo_values(trafo_df, "shift_degree").astype(float) if calculate_voltage_angles else \
359
        np.zeros(len(vnh))
360 1
    if mode == "sc":
361 1
        return vnh, vnl, trafo_shift
362

363 1
    tap_pos = get_trafo_values(trafo_df, "tap_pos")
364 1
    tap_neutral = get_trafo_values(trafo_df, "tap_neutral")
365 1
    tap_diff = tap_pos - tap_neutral
366 1
    tap_phase_shifter = get_trafo_values(trafo_df, "tap_phase_shifter")
367 1
    tap_side = get_trafo_values(trafo_df, "tap_side")
368 1
    tap_step_percent = get_trafo_values(trafo_df, "tap_step_percent")
369 1
    tap_step_degree = get_trafo_values(trafo_df, "tap_step_degree")
370

371 1
    cos = lambda x: np.cos(np.deg2rad(x))
372 1
    sin = lambda x: np.sin(np.deg2rad(x))
373 1
    arctan = lambda x: np.rad2deg(np.arctan(x))
374

375 1
    for side, vn, direction in [("hv", vnh, 1), ("lv", vnl, -1)]:
376 1
        phase_shifters = tap_phase_shifter & (tap_side == side)
377 1
        tap_complex = np.isfinite(tap_step_percent) & np.isfinite(tap_pos) & (tap_side == side) & \
378
                      ~phase_shifters
379 1
        if tap_complex.any():
380 1
            tap_steps = tap_step_percent[tap_complex] * tap_diff[tap_complex] / 100
381 1
            tap_angles = _replace_nan(tap_step_degree[tap_complex])
382 1
            u1 = vn[tap_complex]
383 1
            du = u1 * _replace_nan(tap_steps)
384 1
            vn[tap_complex] = np.sqrt((u1 + du * cos(tap_angles)) ** 2 + (du * sin(tap_angles)) ** 2)
385 1
            trafo_shift[tap_complex] += (arctan(direction * du * sin(tap_angles) /
386
                                                (u1 + du * cos(tap_angles))))
387 1
        if phase_shifters.any():
388 1
            degree_is_set = _replace_nan(tap_step_degree[phase_shifters]) != 0
389 1
            percent_is_set = _replace_nan(tap_step_percent[phase_shifters]) != 0
390 1
            if (degree_is_set & percent_is_set).any():
391 0
                raise UserWarning("Both tap_step_degree and tap_step_percent set for ideal phase shifter")
392 1
            trafo_shift[phase_shifters] += np.where(
393
                (degree_is_set),
394
                (direction * tap_diff[phase_shifters] * tap_step_degree[phase_shifters]),
395
                (direction * 2 * np.rad2deg(np.arcsin(tap_diff[phase_shifters] * \
396
                                                      tap_step_percent[phase_shifters] / 100 / 2)))
397
            )
398 1
    return vnh, vnl, trafo_shift
399

400

401 1
def _replace_nan(array, value=0):
402 1
    mask = np.isnan(array)
403 1
    array[mask] = value
404 1
    return array
405

406 1
def _calc_r_x_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, sn_mva):
407
    """
408
    Calculates (Vectorized) the resitance and reactance according to the
409
    transformer values
410

411
    """
412 1
    parallel = get_trafo_values(trafo_df, "parallel")
413 1
    vk_percent = get_trafo_values(trafo_df, "vk_percent")
414 1
    vkr_percent = get_trafo_values(trafo_df, "vkr_percent")
415 1
    tap_lv = np.square(vn_trafo_lv / vn_lv) * (3* sn_mva)  if mode == 'pf_3ph' else\
416
    np.square(vn_trafo_lv / vn_lv) * sn_mva  # adjust for low voltage side voltage converter
417 1
    sn_trafo_mva = get_trafo_values(trafo_df, "sn_mva")
418 1
    z_sc = vk_percent / 100. / sn_trafo_mva * tap_lv
419 1
    r_sc = vkr_percent / 100. / sn_trafo_mva * tap_lv
420 1
    x_sc = np.sign(z_sc) * np.sqrt(z_sc ** 2 - r_sc ** 2)
421 1
    return r_sc / parallel, x_sc / parallel
422

423

424 1
def _calc_nominal_ratio_from_dataframe(ppc, trafo_df, vn_hv_kv, vn_lv_kv, bus_lookup):
425
    """
426
    Calculates (Vectorized) the off nominal tap ratio::
427

428
                  (vn_hv_kv / vn_lv_kv) / (ub1_in_kv / ub2_in_kv)
429

430
    INPUT:
431
        **net** (Dataframe) - The net for which to calc the tap ratio.
432

433
        **vn_hv_kv** (1d array, float) - The adjusted nominal high voltages
434

435
        **vn_lv_kv** (1d array, float) - The adjusted nominal low voltages
436

437
    OUTPUT:
438
        **tab** (1d array, float) - The off-nominal tap ratio
439
    """
440
    # Calculating tab (trasformer off nominal turns ratio)
441 1
    tap_rat = vn_hv_kv / vn_lv_kv
442 1
    hv_bus = get_trafo_values(trafo_df, "hv_bus")
443 1
    lv_bus = get_trafo_values(trafo_df, "lv_bus")
444 1
    nom_rat = get_values(ppc["bus"][:, BASE_KV], hv_bus, bus_lookup) / \
445
              get_values(ppc["bus"][:, BASE_KV], lv_bus, bus_lookup)
446 1
    return tap_rat / nom_rat
447

448

449 1
def _calc_impedance_parameter(net, ppc):
450 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
451 1
    f, t = net["_pd2ppc_lookups"]["branch"]["impedance"]
452 1
    branch = ppc["branch"]
453 1
    rij, xij, r_asym, x_asym = _calc_impedance_parameters_from_dataframe(net)
454 1
    branch[f:t, BR_R] = rij
455 1
    branch[f:t, BR_X] = xij
456 1
    branch[f:t, BR_R_ASYM] = r_asym
457 1
    branch[f:t, BR_X_ASYM] = x_asym
458 1
    branch[f:t, F_BUS] = bus_lookup[net.impedance["from_bus"].values]
459 1
    branch[f:t, T_BUS] = bus_lookup[net.impedance["to_bus"].values]
460 1
    branch[f:t, BR_STATUS] = net["impedance"]["in_service"].values
461

462

463 1
def _calc_impedance_parameters_from_dataframe(net):
464 1
    impedance = net.impedance
465 1
    sn_impedance = impedance["sn_mva"].values
466 1
    sn_net = net.sn_mva
467 1
    rij = impedance["rft_pu"].values
468 1
    xij = impedance["xft_pu"].values
469 1
    rji = impedance["rtf_pu"].values
470 1
    xji = impedance["xtf_pu"].values
471

472 1
    r = rij / sn_impedance * sn_net
473 1
    x = xij / sn_impedance * sn_net
474 1
    r_asym = (rji - rij) / sn_impedance * sn_net
475 1
    x_asym = (xji - xij) / sn_impedance * sn_net
476 1
    return r, x, r_asym, x_asym
477

478

479 1
def _calc_xward_parameter(net, ppc):
480 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
481 1
    f, t = net["_pd2ppc_lookups"]["branch"]["xward"]
482 1
    branch = ppc["branch"]
483 1
    baseR = np.square(get_values(ppc["bus"][:, BASE_KV], net["xward"]["bus"].values, bus_lookup)) / \
484
            net.sn_mva
485 1
    xw_is = net["_is_elements"]["xward"]
486 1
    branch[f:t, F_BUS] = bus_lookup[net["xward"]["bus"].values]
487 1
    branch[f:t, T_BUS] = bus_lookup[net._pd2ppc_lookups["aux"]["xward"]]
488 1
    branch[f:t, BR_R] = net["xward"]["r_ohm"] / baseR
489 1
    branch[f:t, BR_X] = net["xward"]["x_ohm"] / baseR
490 1
    branch[f:t, BR_STATUS] = xw_is
491

492

493 1
def _gather_branch_switch_info(bus, branch_id, branch_type, net):
494
    # determine at which end the switch is located
495
    # 1 = to-bus/lv-bus; 0 = from-bus/hv-bus
496 1
    branch_id = int(branch_id)
497 1
    lookup = net._pd2ppc_lookups["branch"]
498 1
    if branch_type == "l":
499 1
        side = "to" if net["line"]["to_bus"].at[branch_id] == bus else "from"
500 1
        branch_idx = net["line"].index.get_loc(branch_id)
501 1
        return side, int(bus), int(branch_idx)
502 1
    elif branch_type == "t":
503 1
        side = "hv" if net["trafo"]["hv_bus"].at[branch_id] == bus else "lv"
504 1
        branch_idx = lookup["trafo"][0] + net["trafo"].index.get_loc(branch_id)
505 1
        return side, int(bus), int(branch_idx)
506 1
    elif branch_type == "t3":
507 1
        f, t = lookup["trafo3w"]
508 1
        if net["trafo3w"]["hv_bus"].at[branch_id] == bus:
509 1
            side = "hv"
510 1
            offset = 0
511 1
        elif net["trafo3w"]["mv_bus"].at[branch_id] == bus:
512 1
            side = "mv"
513 1
            offset = (t - f) / 3
514 1
        elif net["trafo3w"]["lv_bus"].at[branch_id] == bus:
515 1
            side = "lv"
516 1
            offset = (t - f) / 3 * 2
517 1
        branch_idx = lookup["trafo3w"][0] + net["trafo3w"].index.get_loc(branch_id) + offset
518 1
        return side, int(bus), int(branch_idx)
519

520

521 1
def _switch_branches(net, ppc):
522 1
    from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX
523 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
524 1
    calculate_voltage_angles = net._options["calculate_voltage_angles"]
525 1
    neglect_open_switch_branches = net._options["neglect_open_switch_branches"]
526 1
    mode = net._options["mode"]
527 1
    open_switches = (net.switch.closed.values == False)
528 1
    n_bus = ppc["bus"].shape[0]
529 1
    for et, element in [("l", "line"), ("t", "trafo"), ("t3", "trafo3w")]:
530 1
        switch_mask = open_switches & (net.switch.et.values == et)
531 1
        if not switch_mask.any():
532 1
            continue
533 1
        nr_open_switches = np.count_nonzero(switch_mask)
534 1
        mapfunc = partial(_gather_branch_switch_info, branch_type=et, net=net)
535 1
        switch_element = net["switch"]["element"].values[switch_mask]
536 1
        switch_buses = net["switch"]["bus"].values[switch_mask]
537 1
        switch_info = np.array(list(map(mapfunc, switch_buses, switch_element)))
538 1
        sw_sides = switch_info[:, 0]
539 1
        sw_bus_index = bus_lookup[switch_info[:, 1].astype(int)]
540 1
        sw_branch_index = switch_info[:, 2].astype(int)
541 1
        if neglect_open_switch_branches:
542
            # deactivate switches which have an open switch instead of creating aux buses
543 0
            ppc["branch"][sw_branch_index, BR_STATUS] = 0
544 0
            continue
545

546 1
        new_buses = np.zeros(shape=(nr_open_switches, ppc["bus"].shape[1]), dtype=float)
547 1
        new_buses[:, :15] = np.array([0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1.1, 0.9, 0, 0])
548 1
        new_indices = np.arange(n_bus, n_bus + nr_open_switches)
549 1
        new_buses[:, 0] = new_indices
550 1
        new_buses[:, BASE_KV] = ppc["bus"][sw_bus_index, BASE_KV]
551 1
        ppc["bus"] = np.vstack([ppc["bus"], new_buses])
552 1
        n_bus += new_buses.shape[0]
553 1
        init_vm = net._options["init_vm_pu"]
554 1
        init_va = net._options["init_va_degree"]
555 1
        for location in np.unique(sw_sides):
556 1
            mask = sw_sides == location
557 1
            buses = new_indices[mask]
558 1
            side = F_BUS if location == "hv" or location == "from" else T_BUS
559 1
            for init, col in [(init_vm, VM), (init_va, VA)]:
560 1
                if isinstance(init, str) and init == "results":
561 1
                    if col == VM:
562 1
                        res_column = net["res_%s" % element]["vm_%s_pu" % location]
563
                    else:
564 1
                        res_column = net["res_%s" % element]["va_%s_degree" % location]
565 1
                    init_values = res_column.loc[switch_element].values[mask]
566
                else:
567 1
                    if element == "line":
568 1
                        opposite_buses = ppc["branch"][sw_branch_index[mask], side].real.astype(int)
569 1
                        init_values = ppc["bus"][opposite_buses, col]
570
                    else:
571 1
                        opposite_side = T_BUS if side == F_BUS else F_BUS
572 1
                        opposite_buses = ppc["branch"][sw_branch_index[mask], opposite_side].real.astype(int)
573 1
                        if col == VM:
574 1
                            taps = ppc["branch"][sw_branch_index[mask], TAP].real
575 1
                            init_values = ppc["bus"][opposite_buses, col] * taps
576
                        else:
577 1
                            if calculate_voltage_angles:
578 1
                                shift = ppc["branch"][sw_branch_index[mask], SHIFT].real.astype(int)
579 1
                                init_values = ppc["bus"][opposite_buses, col] + shift
580
                            else:
581 1
                                init_values = ppc["bus"][opposite_buses, col]
582 1
                ppc["bus"][buses, col] = init_values
583 1
            if mode == "sc":
584 1
                ppc["bus"][buses, C_MAX] = ppc["bus"][opposite_buses, C_MAX]
585 1
                ppc["bus"][buses, C_MIN] = ppc["bus"][opposite_buses, C_MIN]
586 1
            ppc["branch"][sw_branch_index[mask], side] = new_indices[mask]
587

588

589 1
def _branches_with_oos_buses(net, ppc):
590
    """
591
    Updates the ppc["branch"] matrix with the changed from or to values
592
    if the branch is connected to an out of service bus
593

594
    Adds auxiliary buses if branch is connected to an out of service bus
595
    Sets branch out of service if connected to two out of service buses
596

597
    **INPUT**:
598
        **n** - The pandapower format network
599

600
        **ppc** - The PYPOWER format network to fill in values
601
        **bus_is** - The in service buses
602
    """
603 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
604
    # get in service elements
605 1
    _is_elements = net["_is_elements"]
606 1
    bus_is_idx = _is_elements['bus_is_idx']
607 1
    line_is_idx = _is_elements['line_is_idx']
608

609 1
    n_oos_buses = len(net['bus']) - len(bus_is_idx)
610

611
    # only filter lines at oos buses if oos buses exists
612 1
    if n_oos_buses > 0:
613 1
        n_bus = len(ppc["bus"])
614 1
        future_buses = [ppc["bus"]]
615
        # out of service buses
616 1
        bus_oos = np.setdiff1d(net['bus'].index.values, bus_is_idx)
617
        # from buses of line
618 1
        line_buses = net["line"][["from_bus", "to_bus"]].loc[line_is_idx].values
619 1
        f_bus = line_buses[:, 0]
620 1
        t_bus = line_buses[:, 1]
621

622
        # determine on which side of the line the oos bus is located
623 1
        mask_from = np.in1d(f_bus, bus_oos)
624 1
        mask_to = np.in1d(t_bus, bus_oos)
625

626 1
        mask_and = mask_to & mask_from
627 1
        if np.any(mask_and):
628 1
            mask_from[mask_and] = False
629 1
            mask_to[mask_and] = False
630

631
        # get lines that are connected to oos bus at exactly one side
632
        # buses that are connected to two oos buses will be removed by ext2int
633 1
        mask_or = mask_to | mask_from
634
        # check whether buses are connected to line
635 1
        oos_buses_at_lines = np.hstack([f_bus[mask_from], t_bus[mask_to]])
636 1
        n_oos_buses_at_lines = len(oos_buses_at_lines)
637

638
        # only if oos_buses are at lines (they could be isolated as well)
639 1
        if n_oos_buses_at_lines > 0:
640 1
            ls_info = np.zeros((n_oos_buses_at_lines, 3), dtype=int)
641 1
            ls_info[:, 0] = mask_to[mask_or] & ~mask_from[mask_or]
642 1
            ls_info[:, 1] = oos_buses_at_lines
643 1
            ls_info[:, 2] = np.nonzero(np.in1d(net['line'].index, line_is_idx[mask_or]))[0]
644

645
            # ls_info = list(map(mapfunc,
646
            #               line_switches["bus"].values,
647
            #               line_switches["element"].values))
648
            # we now have the following matrix
649
            # 0: 1 if switch is at to_bus, 0 else
650
            # 1: bus of the switch
651
            # 2: position of the line a switch is connected to
652
            # ls_info = np.array(ls_info, dtype=int)
653

654
            # build new buses
655 1
            new_ls_buses = np.zeros(shape=(n_oos_buses_at_lines, ppc["bus"].shape[1]), dtype=float)
656 1
            new_indices = np.arange(n_bus, n_bus + n_oos_buses_at_lines)
657
            # the newly created buses
658 1
            new_ls_buses[:, :15] = np.array([0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1.1, 0.9, 0, 0])
659 1
            new_ls_buses[:, 0] = new_indices
660 1
            new_ls_buses[:, BASE_KV] = get_values(ppc["bus"][:, BASE_KV], ls_info[:, 1], bus_lookup)
661

662 1
            future_buses.append(new_ls_buses)
663

664
            # re-route the end of lines to a new bus
665 1
            ppc["branch"][ls_info[ls_info[:, 0].astype(bool), 2], 1] = \
666
                new_indices[ls_info[:, 0].astype(bool)]
667 1
            ppc["branch"][ls_info[np.logical_not(ls_info[:, 0]), 2], 0] = \
668
                new_indices[np.logical_not(ls_info[:, 0])]
669

670 1
            ppc["bus"] = np.vstack(future_buses)
671

672

673 1
def _calc_switch_parameter(net, ppc):
674
    """
675
    calculates the line parameter in per unit.
676

677
    **INPUT**:
678
        **net** -The pandapower format network
679

680
    **RETURN**:
681
        **t** - Temporary line parameter. Which is a complex128
682
                Nunmpy array. with the following order:
683
                0:bus_a; 1:bus_b; 2:r_pu; 3:x_pu; 4:b_pu
684
    """
685 1
    rx_ratio = net["_options"]["switch_rx_ratio"]
686 1
    rz_ratio = rx_ratio / np.sqrt(1 + rx_ratio ** 2)
687 1
    xz_ratio = 1 / np.sqrt(1 + rx_ratio ** 2)
688

689 1
    f, t = net["_pd2ppc_lookups"]["branch"]["switch"]
690 1
    branch = ppc["branch"]
691 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
692 1
    switch = net.switch[net._impedance_bb_switches]
693 1
    fb = bus_lookup[switch["bus"].values]
694 1
    tb = bus_lookup[switch["element"].values]
695 1
    baseR = np.square(ppc["bus"][fb, BASE_KV]) / net.sn_mva
696 1
    branch[f:t, F_BUS] = fb
697 1
    branch[f:t, T_BUS] = tb
698

699 1
    z_switch = switch['z_ohm'].values
700
    # x_switch will have the same value of r_switch to avoid zero dividence
701 1
    branch[f:t, BR_R] = z_switch / baseR * rz_ratio
702 1
    branch[f:t, BR_X] = z_switch / baseR * xz_ratio
703

704

705 1
def _end_temperature_correction_factor(net, short_circuit=False):
706
    """
707
    Function to calculate resistance correction factor for the given temperature ("endtemp_degree").
708
    When multiplied by the factor, the value of r_ohm_per_km will correspond to the resistance at
709
    the given temperature.
710

711
    In case of short circuit calculation, the relevant value for the temperature is
712
    "endtemp_degree", which stands for the final temperature of a line after the short circuit.
713
    The temperature coefficient "alpha" is a constant value of 0.004 in the short circuit
714
    calculation standard IEC 60909-0:2016.
715

716
    In case of a load flow calculation, the relelvant parameter is "temperature_degree_celsius",
717
    which is specified by the user and allows calculating load flow for a given operating
718
    temperature.
719

720
    The alpha value can be provided according to the used material for the purpose of load flow
721
    calculation, e.g. 0.0039 for copper or 0.00403 for aluminum. If alpha is not provided in the
722
    net.line table, the default value of 0.004 is used.
723

724
    The calculation of the electrical resistance is based on the formula R = R20(1+alpha*(T-20°C)),
725
    where R is the calculated resistance, R20 is the resistance at 20 °C, alpha is the temperature
726
    coefficient of resistance of the conducting material and T is the line temperature in °C.
727
    Accordingly, the resulting correction factor is (1+alpha*(T-20°C)).
728

729
    Args:
730
        net: pandapowerNet
731
        short_circuit: whether the factor is calculated in the scope of a short circuit calculation
732

733
    Returns:
734
        correction factor for line R, by which the line parameter should be multiplied to
735
                obtain the value of resistance at line temperature "endtemp_degree"
736

737
    """
738

739 1
    if short_circuit:
740
        # endtemp_degree is line temperature that is reached as the result of a short circuit
741
        # this value is the property of the lines
742 1
        if "endtemp_degree" not in net.line.columns:
743 0
            raise UserWarning("Specify end temperature for lines in net.line.endtemp_degree")
744

745 1
        delta_t_degree_celsius = net.line.endtemp_degree.values.astype(np.float64) - 20
746
        # alpha is the temperature correction factor for the electric resistance of the material
747
        # formula from standard, used this way in short-circuit calculation
748 1
        alpha = 4e-3
749
    else:
750
        # temperature_degree_celsius is line temperature for load flow calculation
751 1
        if "temperature_degree_celsius" not in net.line.columns:
752 0
            raise UserWarning("Specify line temperature in net.line.temperature_degree_celsius")
753

754 1
        delta_t_degree_celsius = net.line.temperature_degree_celsius.values.astype(np.float64) - 20
755

756 1
        if 'alpha' in net.line.columns:
757 1
            alpha = net.line.alpha.values.astype(np.float64)
758
        else:
759 1
            alpha = 4e-3
760

761 1
    r_correction_for_temperature = 1 + alpha * delta_t_degree_celsius
762

763 1
    return r_correction_for_temperature
764

765

766 1
def _transformer_correction_factor(vk, vkr, sn, cmax):
767
    """
768
        2W-Transformer impedance correction factor in short circuit calculations,
769
        based on the IEC 60909-0:2016 standard.
770
        Args:
771
            vk: transformer short-circuit voltage, percent
772
            vkr: real-part of transformer short-circuit voltage, percent
773
            sn: transformer rating, kVA
774
            cmax: voltage factor to account for maximum worst-case currents, based on the lv side
775

776
        Returns:
777
            kt: transformer impedance correction factor for short-circuit calculations
778
        """
779 1
    zt = vk / 100 / sn
780 1
    rt = vkr / 100 / sn
781 1
    xt = np.sqrt(zt ** 2 - rt ** 2)
782 1
    kt = 0.95 * cmax / (1 + .6 * xt * sn)
783 1
    return kt
784

785

786 1
def get_is_lines(net):
787
    """
788
    get indices of lines that are in service and save that information in net
789
    """
790 0
    _is_elements = net["_is_elements"]
791 0
    _is_elements["line"] = net["line"][net["line"]["in_service"].values.astype(bool)]
792

793

794 1
def _trafo_df_from_trafo3w(net):
795 1
    nr_trafos = len(net["trafo3w"])
796 1
    trafo2 = dict()
797 1
    sides = ["hv", "mv", "lv"]
798 1
    mode = net._options["mode"]
799 1
    loss_side = net._options["trafo3w_losses"].lower()
800 1
    nr_trafos = len(net["trafo3w"])
801 1
    t3 = net["trafo3w"]
802 1
    _calculate_sc_voltages_of_equivalent_transformers(t3, trafo2, mode)
803 1
    _calculate_3w_tap_changers(t3, trafo2, sides)
804 1
    zeros = np.zeros(len(net.trafo3w))
805 1
    aux_buses = net._pd2ppc_lookups["aux"]["trafo3w"]
806 1
    trafo2["hv_bus"] = {"hv": t3.hv_bus.values, "mv": aux_buses, "lv": aux_buses}
807 1
    trafo2["lv_bus"] = {"hv": aux_buses, "mv": t3.mv_bus.values, "lv": t3.lv_bus.values}
808 1
    trafo2["in_service"] = {side: t3.in_service.values for side in sides}
809 1
    trafo2["i0_percent"] = {side: t3.i0_percent.values if loss_side == side else zeros for side in sides}
810 1
    trafo2["pfe_kw"] = {side: t3.pfe_kw.values if loss_side == side else zeros for side in sides}
811 1
    trafo2["vn_hv_kv"] = {side: t3.vn_hv_kv.values for side in sides}
812 1
    trafo2["vn_lv_kv"] = {side: t3["vn_%s_kv" % side].values for side in sides}
813 1
    trafo2["shift_degree"] = {"hv": np.zeros(nr_trafos), "mv": t3.shift_mv_degree.values,
814
                              "lv": t3.shift_lv_degree.values}
815 1
    trafo2["tap_phase_shifter"] = {side: np.zeros(nr_trafos).astype(bool) for side in sides}
816 1
    trafo2["parallel"] = {side: np.ones(nr_trafos) for side in sides}
817 1
    trafo2["df"] = {side: np.ones(nr_trafos) for side in sides}
818 1
    if net._options["mode"] == "opf" and "max_loading_percent" in net.trafo3w:
819 1
        trafo2["max_loading_percent"] = {side: net.trafo3w.max_loading_percent.values for side in sides}
820 1
    return {var: np.concatenate([trafo2[var][side] for side in sides]) for var in trafo2.keys()}
821

822

823 1
def _calculate_sc_voltages_of_equivalent_transformers(t3, t2, mode):
824 1
    vk_3w = np.stack([t3.vk_hv_percent.values, t3.vk_mv_percent.values, t3.vk_lv_percent.values])
825 1
    vkr_3w = np.stack([t3.vkr_hv_percent.values, t3.vkr_mv_percent.values, t3.vkr_lv_percent.values])
826 1
    sn = np.stack([t3.sn_hv_mva.values, t3.sn_mv_mva.values, t3.sn_lv_mva.values])
827

828 1
    vk_2w_delta = z_br_to_bus_vector(vk_3w, sn)
829 1
    vkr_2w_delta = z_br_to_bus_vector(vkr_3w, sn)
830 1
    if mode == "sc":
831 1
        kt = _transformer_correction_factor(vk_3w, vkr_3w, sn, 1.1)
832 1
        vk_2w_delta *= kt
833 1
        vkr_2w_delta *= kt
834 1
    vki_2w_delta = np.sqrt(vk_2w_delta ** 2 - vkr_2w_delta ** 2)
835 1
    vkr_2w = wye_delta_vector(vkr_2w_delta, sn)
836 1
    vki_2w = wye_delta_vector(vki_2w_delta, sn)
837 1
    vk_2w = np.sign(vki_2w) * np.sqrt(vki_2w ** 2 + vkr_2w ** 2)
838 1
    if np.any(vk_2w == 0):
839 0
        raise UserWarning("Equivalent transformer with zero impedance!")
840 1
    t2["vk_percent"] = {"hv": vk_2w[0, :], "mv": vk_2w[1, :], "lv": vk_2w[2, :]}
841 1
    t2["vkr_percent"] = {"hv": vkr_2w[0, :], "mv": vkr_2w[1, :], "lv": vkr_2w[2, :]}
842 1
    t2["sn_mva"] = {"hv": sn[0, :], "mv": sn[1, :], "lv": sn[2, :]}
843

844

845 1
def z_br_to_bus_vector(z, sn):
846 1
    return sn[0, :] * np.array([z[0, :] / sn[[0, 1], :].min(axis=0), z[1, :] /
847
                                sn[[1, 2], :].min(axis=0), z[2, :] / sn[[0, 2], :].min(axis=0)])
848

849

850 1
def wye_delta(zbr_n, s):
851 0
    return .5 * s / s[0] * np.array([(zbr_n[0] + zbr_n[2] - zbr_n[1]),
852
                                     (zbr_n[1] + zbr_n[0] - zbr_n[2]),
853
                                     (zbr_n[2] + zbr_n[1] - zbr_n[0])])
854

855

856 1
def wye_delta_vector(zbr_n, s):
857 1
    return .5 * s / s[0, :] * np.array([(zbr_n[0, :] + zbr_n[2, :] - zbr_n[1, :]),
858
                                        (zbr_n[1, :] + zbr_n[0, :] - zbr_n[2, :]),
859
                                        (zbr_n[2, :] + zbr_n[1, :] - zbr_n[0, :])])
860

861

862 1
def _calculate_3w_tap_changers(t3, t2, sides):
863 1
    tap_variables = ["tap_side", "tap_pos", "tap_neutral", "tap_max", "tap_min", "tap_step_percent",
864
                     "tap_step_degree"]
865 1
    sides = ["hv", "mv", "lv"]
866 1
    nr_trafos = len(t3)
867 1
    empty = np.zeros(nr_trafos)
868 1
    empty.fill(np.nan)
869 1
    tap_arrays = {var: {side: empty.copy() for side in sides} for var in tap_variables}
870 1
    tap_arrays["tap_side"] = {side: np.array([None] * nr_trafos) for side in sides}
871 1
    at_star_point = t3.tap_at_star_point.values
872 1
    any_at_star_point = at_star_point.any()
873 1
    for side in sides:
874 1
        tap_mask = t3.tap_side.values == side
875 1
        for var in tap_variables:
876 1
            tap_arrays[var][side][tap_mask] = t3[var].values[tap_mask]
877

878
        # t3 trafos with tap changer at terminals
879 1
        tap_arrays["tap_side"][side][tap_mask] = "hv" if side == "hv" else "lv"
880

881
        # t3 trafos with tap changer at star points
882 1
        if any_at_star_point:
883 0
            mask_star_point = tap_mask & at_star_point
884 0
            tap_arrays["tap_side"][side][mask_star_point] = "lv" if side == "hv" else "hv"
885 0
            tap_arrays["tap_step_degree"][side][mask_star_point] += 180
886 1
    t2.update(tap_arrays)

Read our documentation on viewing source code .

Loading