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
import numpy as np
8 1
import pandas as pd
9

10 1
from pandapower.auxiliary import _sum_by_group
11 1
from pandapower.pypower.idx_bus import BASE_KV
12 1
from pandapower.pypower.idx_gen import GEN_BUS, MBASE
13 1
from pandapower.shortcircuit.idx_brch import IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T
14 1
from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX, KAPPA, R_EQUIV, IKSS1, IP, ITH, X_EQUIV, IKSS2, IKCV, M
15 1
from pandapower.shortcircuit.impedance import _calc_zbus_diag
16

17 1
from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
18 1
from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
19

20 1
def _calc_ikss(net, ppc, bus=None):
21
    # Vectorized for multiple bus
22 1
    if bus is None:
23
        # Slice(None) is equal to : select
24 1
        bus_idx = slice(None)
25
    else:
26 1
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
27

28 1
    fault = net._options["fault"]
29 1
    case = net._options["case"]
30 1
    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
31 1
    ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
32

33 1
    z_equiv = abs(ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j)
34 1
    if fault == "3ph":
35 1
        ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
36 1
    elif fault == "2ph":
37 1
        ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / 2 * ppc["baseMVA"]
38 1
    _current_source_current(net, ppc)
39

40

41 1
def _calc_ikss_1ph(net, ppc, ppc_0, bus=None):
42
    # Vectorized for multiple bus
43 1
    if bus is None:
44
        # Slice(None) is equal to : select
45 1
        bus_idx = slice(None)
46
    else:
47 1
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
48

49 1
    case = net._options["case"]
50 1
    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
51 1
    ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
52 1
    ppc_0["internal"]["baseI"] = ppc_0["bus"][:, BASE_KV] * np.sqrt(3) / ppc_0["baseMVA"]
53

54 1
    z_equiv = abs((ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j) * 2 +
55
                  (ppc_0["bus"][bus_idx, R_EQUIV] + ppc_0["bus"][bus_idx, X_EQUIV] * 1j))
56

57 1
    ppc_0["bus"][bus_idx, IKSS1] = c / z_equiv / ppc_0["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
58 1
    ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc["baseMVA"]
59 1
    _current_source_current(net, ppc)
60

61

62 1
def _current_source_current(net, ppc):
63 1
    ppc["bus"][:, IKCV] = 0
64 1
    ppc["bus"][:, IKSS2] = 0
65 1
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
66 1
    if not False in net.sgen.current_source.values:
67 1
        sgen = net.sgen[net._is_elements["sgen"]]
68
    else:
69 0
        sgen = net.sgen[net._is_elements["sgen"] & net.sgen.current_source]
70 1
    if len(sgen) == 0:
71 1
        return
72 1
    if any(pd.isnull(sgen.sn_mva)):
73 0
        raise ValueError("sn_mva needs to be specified for all sgens in net.sgen.sn_mva")
74 1
    baseI = ppc["internal"]["baseI"]
75 1
    sgen_buses = sgen.bus.values
76 1
    sgen_buses_ppc = bus_lookup[sgen_buses]
77 1
    if not "k" in sgen:
78 0
        raise ValueError("Nominal to short-circuit current has to specified in net.sgen.k")
79 1
    i_sgen_pu = sgen.sn_mva.values / net.sn_mva * sgen.k.values
80 1
    buses, ikcv_pu, _ = _sum_by_group(sgen_buses_ppc, i_sgen_pu, i_sgen_pu)
81 1
    ppc["bus"][buses, IKCV] = ikcv_pu
82 1
    if net["_options"]["inverse_y"]:
83 1
        Zbus = ppc["internal"]["Zbus"]
84 1
        ppc["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppc["bus"][:, IKCV] * -1j) / baseI)
85
    else:
86 1
        ybus_fact = ppc["internal"]["ybus_fact"]
87 1
        diagZ = _calc_zbus_diag(net, ppc)
88 1
        ppc["bus"][:, IKSS2] = abs(ybus_fact(ppc["bus"][:, IKCV] * -1j) / diagZ / baseI)
89 1
    ppc["bus"][buses, IKCV] /= baseI[buses]
90

91

92 1
def _calc_ip(net, ppc):
93 1
    ip = np.sqrt(2) * (ppc["bus"][:, KAPPA] * ppc["bus"][:, IKSS1] + ppc["bus"][:, IKSS2])
94 1
    ppc["bus"][:, IP] = ip
95

96

97 1
def _calc_ith(net, ppc):
98 1
    tk_s = net["_options"]["tk_s"]
99 1
    kappa = ppc["bus"][:, KAPPA]
100 1
    f = 50
101 1
    n = 1
102 1
    m = (np.exp(4 * f * tk_s * np.log(kappa - 1)) - 1) / (2 * f * tk_s * np.log(kappa - 1))
103 1
    m[np.where(kappa > 1.99)] = 0
104 1
    ppc["bus"][:, M] = m
105 1
    ith = (ppc["bus"][:, IKSS1] + ppc["bus"][:, IKSS2]) * np.sqrt(m + n)
106 1
    ppc["bus"][:, ITH] = ith
107

108

109 1
def _calc_ib_generator(net, ppci):
110 0
    Zbus = ppci["internal"]["Zbus"]
111 0
    baseI = ppci["internal"]["baseI"]
112 0
    tk_s = net._options['tk_s']
113 0
    c = 1.1
114

115 0
    z_equiv = ppci["bus"][:, R_EQUIV] + ppci["bus"][:, X_EQUIV] * 1j
116 0
    I_ikss = c / z_equiv / ppci["bus"][:, BASE_KV] / np.sqrt(3) * ppci["baseMVA"]
117

118
    # calculate voltage source branch current
119
    # I_ikss = ppci["bus"][:, IKSS1]
120 0
    V_ikss = (I_ikss * baseI) * Zbus
121

122 0
    gen = net["gen"][net._is_elements["gen"]]
123 0
    gen_vn_kv = gen.vn_kv.values
124

125 0
    gen_buses = ppci['gen'][:, GEN_BUS].astype(np.int64)
126 0
    gen_mbase = ppci['gen'][:, MBASE]
127 0
    gen_i_rg = gen_mbase / (np.sqrt(3) * gen_vn_kv)
128

129 0
    gen_buses_ppc, gen_sn_mva, I_rG = _sum_by_group(gen_buses, gen_mbase, gen_i_rg)
130

131
    # shunt admittance of generator buses and generator short circuit current
132
    # YS = ppci["bus"][gen_buses_ppc, GS] + ppci["bus"][gen_buses_ppc, BS] * 1j
133
    # I_kG = V_ikss.T[:, gen_buses_ppc] * YS / baseI[gen_buses_ppc]
134

135 0
    xdss_pu = gen.xdss_pu.values
136 0
    rdss_pu = gen.rdss_pu.values
137 0
    cosphi = gen.cos_phi.values
138 0
    X_dsss = xdss_pu * np.square(gen_vn_kv) / gen_mbase
139 0
    R_dsss = rdss_pu * np.square(gen_vn_kv) / gen_mbase
140

141 0
    K_G = ppci['bus'][gen_buses, BASE_KV] / gen_vn_kv * c / (1 + xdss_pu * np.sin(np.arccos(cosphi)))
142 0
    Z_G = (R_dsss + 1j * X_dsss)
143

144 0
    I_kG = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3) / (Z_G * K_G) * ppci["baseMVA"]
145

146 0
    dV_G = 1j * X_dsss * K_G * I_kG
147 0
    V_Is = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3)
148

149
    # I_kG_contribution = I_kG.sum(axis=1)
150
    # ratio_SG_ikss = I_kG_contribution / I_ikss
151
    # close_to_SG = ratio_SG_ikss > 5e-2
152

153 0
    close_to_SG = I_kG / I_rG > 2
154

155 0
    if tk_s == 2e-2:
156 0
        mu = 0.84 + 0.26 * np.exp(-0.26 * abs(I_kG) / I_rG)
157 0
    elif tk_s == 5e-2:
158 0
        mu = 0.71 + 0.51 * np.exp(-0.3 * abs(I_kG) / I_rG)
159 0
    elif tk_s == 10e-2:
160 0
        mu = 0.62 + 0.72 * np.exp(-0.32 * abs(I_kG) / I_rG)
161 0
    elif tk_s >= 25e-2:
162 0
        mu = 0.56 + 0.94 * np.exp(-0.38 * abs(I_kG) / I_rG)
163
    else:
164 0
        raise UserWarning('not implemented for other tk_s than 20ms, 50ms, 100ms and >=250ms')
165

166 0
    mu = np.clip(mu, 0, 1)
167

168 0
    I_ikss_G = abs(I_ikss - np.sum((1 - mu) * I_kG, axis=1))
169

170
    # I_ikss_G = I_ikss - np.sum(abs(V_ikss.T[:, gen_buses_ppc]) * (1-mu) * I_kG, axis=1)
171

172 0
    I_ikss_G = abs(I_ikss - np.sum(dV_G / V_Is * (1 - mu) * I_kG, axis=1))
173

174 0
    return I_ikss_G
175

176

177 1
def _calc_single_bus_sc(net, ppc, bus):
178
    # case = net._options["case"]
179 1
    bus_idx = net._pd2ppc_lookups["bus"][bus]
180 1
    n = ppc["bus"].shape[0]
181 1
    Zbus = ppc["internal"]["Zbus"]
182
    #    Yf = ppc["internal"]["Yf"]
183
    #    Yt = ppc["internal"]["Yf"]
184 1
    baseI = ppc["internal"]["baseI"]
185
    #    fb = np.real(ppc["branch"][:, 0]).astype(int)
186
    #    tb = np.real(ppc["branch"][:, 1]).astype(int)
187
    # c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
188

189
    # calculate voltage source branch current
190 1
    V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
191 1
    V = V_ikss[:, bus_idx]
192
    #    ikss_all_f = np.conj(Yf.dot(V_ikss))
193
    #    ikss_all_t = np.conj(Yt.dot(V_ikss))
194 1
    current_sources = any(ppc["bus"][:, IKCV]) > 0
195 1
    if current_sources:
196 1
        current = np.tile(-ppc["bus"][:, IKCV], (n, 1))
197 1
        np.fill_diagonal(current, current.diagonal() + ppc["bus"][:, IKSS2])
198 1
        V_source = np.dot((current * baseI), Zbus).T
199 1
        V = V + V_source[:, bus_idx]
200
    # add current source branch current if there is one
201
    #    ppc["branch"][:, IKSS_F] = abs(ikss_all_f[:, bus_idx] / baseI[fb])
202
    #    ppc["branch"][:, IKSS_T] = abs(ikss_all_t[:, bus_idx] / baseI[tb])
203 1
    calc_branch_results(net, ppc, V)
204

205 1
def _calc_single_bus_sc_no_y_inv(net, ppc, bus):
206
    # Vectorized for multiple bus
207 1
    if bus is None:
208
        # Slice(None) is equal to : select
209 0
        bus_idx = slice(None)
210
    else:
211 1
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
212

213 1
    ybus = ppc["internal"]["Ybus"]
214 1
    ybus_fact = ppc["internal"]["ybus_fact"]
215
    # case = net._options["case"]
216 1
    baseI = ppc["internal"]["baseI"]
217
    # vqj = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX] #this is the source voltage in per unit (VQj)
218

219
    # Solve Ikss from voltage source
220 1
    n_bus = ybus.shape[0]
221

222
    # ybus_sub_mask = (np.arange(ybus.shape[0]) != bus_idx)
223
    # V_ikss = np.zeros(n_bus, dtype=np.complex)
224
    # V_ikss[bus_idx] = vqj[bus_idx]
225

226
    # # Solve Ax = b
227
    # b = np.zeros(n_bus-1, dtype=np.complex) -\
228
    #     (ybus[:, ~ybus_sub_mask].toarray())[ybus_sub_mask].ravel() * V_ikss[bus_idx]
229
    # ybus_sub = ybus[ybus_sub_mask, :][:, ybus_sub_mask]
230
    # x = spsolve(ybus_sub, b)
231

232
    # V_ikss[ybus_sub_mask] = x
233
    # I_ikss = np.zeros(n_bus, dtype=np.complex)
234
    # I_ikss[bus_idx] = np.dot(ybus[bus_idx, :].toarray(), V_ikss)
235
    # V = V_ikss
236

237
    # Version 2
238 1
    I_ikss = np.zeros(n_bus, dtype=np.complex)
239 1
    I_ikss[bus_idx] = ppc["bus"][bus_idx, IKSS1]
240 1
    V_ikss = ybus_fact(I_ikss * baseI)
241 1
    V = V_ikss
242

243
    #TODO include current sources
244 1
    current_sources = any(ppc["bus"][:, IKCV]) > 0
245 1
    if current_sources:
246 1
        current = -ppc["bus"][:, IKCV]
247 1
        current[bus_idx] += ppc["bus"][bus_idx, IKSS2]
248 1
        V_source = ybus_fact(current)
249 1
        V += V_source
250

251 1
    calc_branch_results(net, ppc, V)
252

253

254 1
def calc_branch_results(net, ppci, V):
255 1
    Ybus = ppci["internal"]["Ybus"]
256 1
    Yf = ppci["internal"]["Yf"]
257 1
    Yt = ppci["internal"]["Yt"]
258 1
    baseMVA, bus, gen, branch, ref, _, _, _, _, _, ref_gens = _get_pf_variables_from_ppci(ppci)
259 1
    bus, gen, branch = pfsoln_pypower(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens)
260 1
    ppci["bus"], ppci["gen"], ppci["branch"] = bus, gen, branch
261

262

263 1
def _calc_branch_currents(net, ppc, bus):
264
    # Vectorized for multiple bus
265 1
    if bus is None:
266
        # Slice(None) is equal to select all
267 1
        bus = net.bus.index
268

269 1
    bus_idx = net._pd2ppc_lookups["bus"][bus]
270
    # Select only in service bus for sc calculation
271 1
    bus_idx = bus_idx[bus_idx < ppc['bus'].shape[0]]
272 1
    n_sc_bus = np.shape(bus_idx)[0]
273

274 1
    case = net._options["case"]
275

276 1
    Yf = ppc["internal"]["Yf"]
277 1
    Yt = ppc["internal"]["Yt"]
278 1
    baseI = ppc["internal"]["baseI"]
279 1
    n_bus = ppc["bus"].shape[0]
280 1
    fb = np.real(ppc["branch"][:, 0]).astype(int)
281 1
    tb = np.real(ppc["branch"][:, 1]).astype(int)
282 1
    minmax = np.nanmin if case == "min" else np.nanmax
283

284
    # calculate voltage source branch current
285 1
    if net["_options"]["inverse_y"]:
286 1
        Zbus = ppc["internal"]["Zbus"]
287 1
        V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
288 1
        V_ikss = V_ikss[:, bus_idx]
289
    else:
290 1
        ybus_fact = ppc["internal"]["ybus_fact"]
291 1
        V_ikss = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
292 1
        for ix, b in enumerate(bus_idx):
293 1
            ikss = np.zeros(n_bus, dtype=np.complex)
294 1
            ikss[b] = ppc["bus"][b, IKSS1] * baseI[b]
295 1
            V_ikss[:, ix] = ybus_fact(ikss)
296

297 1
    ikss1_all_f = np.conj(Yf.dot(V_ikss))
298 1
    ikss1_all_t = np.conj(Yt.dot(V_ikss))
299 1
    ikss1_all_f[abs(ikss1_all_f) < 1e-10] = 0.
300 1
    ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.
301

302
    # add current source branch current if there is one
303 1
    current_sources = any(ppc["bus"][:, IKCV]) > 0
304 1
    if current_sources:
305 1
        current = np.tile(-ppc["bus"][:, IKCV], (n_sc_bus, 1))
306 1
        for ix, b in enumerate(bus_idx):
307 1
            current[ix, b] += ppc["bus"][b, IKSS2]
308

309
        # calculate voltage source branch current
310 1
        if net["_options"]["inverse_y"]:
311 1
            Zbus = ppc["internal"]["Zbus"]
312 1
            V = np.dot((current * baseI), Zbus).T
313
        else:
314 1
            ybus_fact = ppc["internal"]["ybus_fact"]
315 1
            V = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
316 1
            for ix, b in enumerate(bus_idx):
317 1
                V[:, ix] = ybus_fact(current[ix, :] * baseI[b])
318

319 1
        fb = np.real(ppc["branch"][:, 0]).astype(int)
320 1
        tb = np.real(ppc["branch"][:, 1]).astype(int)
321 1
        ikss2_all_f = np.conj(Yf.dot(V))
322 1
        ikss2_all_t = np.conj(Yt.dot(V))
323

324 1
        ikss_all_f = abs(ikss1_all_f + ikss2_all_f)
325 1
        ikss_all_t = abs(ikss1_all_t + ikss2_all_t)
326
    else:
327 1
        ikss_all_f = abs(ikss1_all_f)
328 1
        ikss_all_t = abs(ikss1_all_t)
329

330 1
    if net._options["return_all_currents"]:
331 1
        ppc["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
332 1
        ppc["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
333
    else:
334 1
        ikss_all_f[abs(ikss_all_f) < 1e-10] = np.nan
335 1
        ikss_all_t[abs(ikss_all_t) < 1e-10] = np.nan
336 1
        ppc["branch"][:, IKSS_F] = np.nan_to_num(minmax(ikss_all_f, axis=1) / baseI[fb])
337 1
        ppc["branch"][:, IKSS_T] = np.nan_to_num(minmax(ikss_all_t, axis=1) / baseI[tb])
338

339 1
    if net._options["ip"]:
340 1
        kappa = ppc["bus"][:, KAPPA]
341 1
        if current_sources:
342 1
            ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa[bus_idx] + ikss2_all_f)
343 1
            ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa[bus_idx] + ikss2_all_t)
344
        else:
345 1
            ip_all_f = np.sqrt(2) * ikss1_all_f * kappa[bus_idx]
346 1
            ip_all_t = np.sqrt(2) * ikss1_all_t * kappa[bus_idx]
347

348 1
        if net._options["return_all_currents"]:
349 1
            ppc["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
350 1
            ppc["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
351
        else:
352 1
            ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
353 1
            ip_all_t[abs(ip_all_t) < 1e-10] = np.nan
354 1
            ppc["branch"][:, IP_F] = np.nan_to_num(minmax(abs(ip_all_f), axis=1) / baseI[fb])
355 1
            ppc["branch"][:, IP_T] = np.nan_to_num(minmax(abs(ip_all_t), axis=1) / baseI[tb])
356

357 1
    if net._options["ith"]:
358 1
        n = 1
359 1
        m = ppc["bus"][bus_idx, M]
360 1
        ith_all_f = ikss_all_f * np.sqrt(m + n)
361 1
        ith_all_t = ikss_all_t * np.sqrt(m + n)
362

363 1
        if net._options["return_all_currents"]:
364 1
            ppc["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
365 1
            ppc["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
366
        else:
367 1
            ppc["branch"][:, ITH_F] = np.nan_to_num(minmax(ith_all_f, axis=1) / baseI[fb])
368 1
            ppc["branch"][:, ITH_T] = np.nan_to_num(minmax(ith_all_t, axis=1) / baseI[fb])

Read our documentation on viewing source code .

Loading