IEC standard specifies usage of rated voltage of the node to do short circuit calculations. This is because it makes no assumption on actual voltage but assumes that nominally it would be 1.0 p.u. for a bus and at extreme, 1.1.
PGM is not 100% conformant with IEC 60909.
We do not calculate short circuit currents using rated voltage of fault nodes: c * u_rated. But we use c * u_rated * k_eq which includes any transformation ratio k_eq involved from source to the fault node. This k factor is statically present in the grid and not dependent on load conditions hence it is present in exact form in the short circuit situation too.
The design decision was that since our approach would be more accurate compared to the standard, we do not have to stick with certain assumptions.
The other decisions like inclusion of magnetisation current, tap positions and capcitance in the short circuit calculations are easily managed on the user side by setting i0, c1, c0, tap_pos, tap_nom as 0. In the case of voltage, determining k_eq especially with multiple sources becomes more complex to bypass.
Many of MV grids have varied transformer secondary ratings compared to a fixed 10kV node voltage rating, hence this issue concerns those cases.
As expected, we also deviate from other calculation engines in this regard (See in the script below). Hence our validation does not become straightforward.
Given these points, let us revisit this decision again.
Comparison script
Reproduces PGM calculations, calculations inside PGM replicated by hand, hand calculations via IEC standard, vision (Calculated externally) and pandapower. (TLDR to relevant functions or output directly.)
The parameters are unrealistic but suitable for understanding the calculation.
source -- node_1 -- transformer -- node_2
import warnings
from pprint import pprint
import numpy as np
from power_grid_model import (
ComponentType,
CalculationMethod,
PowerGridModel,
initialize_array,
DatasetType,
AttributeType as AT,
FaultType,
ShortCircuitVoltageScaling,
FaultPhase
)
import pandapower as pp
# Constants
s_base = 1e6
u_rated_source_node = 100e3
z_base_source_node = u_rated_source_node * u_rated_source_node / s_base
i_base_source_node = s_base / u_rated_source_node
# u_rated_transformer_node will be set in the function calls below. Not a constant
# Transformer params
u1 = 100e3
u2 = 10e3
sn = 1000e6
uk = 0.1
def calculate_pgm_engine(skss, c, fault_node, u_rated_transformer_node):
skss_adj = skss / c # Assuming #1352 is resolved
# node
node = initialize_array(DatasetType.input, ComponentType.node, 2)
node[AT.id] = [1, 2]
node[AT.u_rated] = [u_rated_source_node, u_rated_transformer_node]
# source
source = initialize_array(DatasetType.input, ComponentType.source, 1)
source[AT.id] = [3]
source[AT.node] = [1]
source[AT.status] = [1]
source[AT.u_ref] = [1.0]
source[AT.sk] = [skss_adj]
source[AT.rx_ratio] = [1e-4]
# transformer
transformer = initialize_array(DatasetType.input, ComponentType.transformer, 1)
transformer[AT.id] = [4]
transformer[AT.from_node] = [1]
transformer[AT.to_node] = [2]
transformer[AT.from_status] = [1]
transformer[AT.to_status] = [1]
transformer[AT.u1] = [u1]
transformer[AT.u2] = [u2]
transformer[AT.sn] = [sn]
transformer[AT.uk] = [uk]
transformer[AT.pk] = [0.0]
transformer[AT.p0] = [0.0]
transformer[AT.i0] = [0.0]
transformer[AT.tap_side] = [0]
transformer[AT.tap_pos] = [0]
transformer[AT.tap_nom] = [0]
transformer[AT.tap_min] = [-10]
transformer[AT.tap_max] = [10]
transformer[AT.tap_size] = [1000]
transformer[AT.clock] = [0]
transformer[AT.winding_from] = [0]
transformer[AT.winding_to] = [0]
# fault
fault = initialize_array(DatasetType.input, ComponentType.fault, 1)
fault[AT.id] = [5]
fault[AT.fault_object] = [fault_node]
fault[AT.status] = [1]
fault[AT.fault_type] = FaultType.three_phase
fault[AT.fault_phase] = FaultPhase.abc
fault[AT.r_f] = [0.0]
fault[AT.x_f] = [0.0]
# all
input_data = {
ComponentType.node: node,
ComponentType.source: source,
ComponentType.transformer: transformer,
ComponentType.fault: fault,
}
model = PowerGridModel(input_data)
short_circuit_voltage_scaling = (
ShortCircuitVoltageScaling.minimum
if np.isclose(c, 1.0)
else ShortCircuitVoltageScaling.maximum
)
output_data = model.calculate_short_circuit(
calculation_method=CalculationMethod.iec60909,
short_circuit_voltage_scaling=short_circuit_voltage_scaling,
)
return {
"i_f": output_data[ComponentType.fault][AT.i_f][0, 0]
}
def calculate_pandapower(skss, c, fault_node, u_rated_transformer_node):
with warnings.catch_warnings():
warnings.simplefilter("ignore")
#create empty net
net = pp.create_empty_network()
#create buses
bus1 = pp.create_bus(net, vn_kv=u_rated_source_node, name="Bus 1")
bus2 = pp.create_bus(net, vn_kv=u_rated_transformer_node, name="Bus 2")
pp.create_ext_grid(net, bus=bus1, vm_pu=1.0, s_sc_max_mva=skss, s_sc_min_mva=skss, rx_max=1e-4, rx_min=1e-4)
pp.create_transformer_from_parameters(net, hv_bus=bus1, lv_bus=bus2, sn_mva=sn, vn_hv_kv=u1, vn_lv_kv=u2, vk_percent=uk*100, vkr_percent=0.0, pfe_kw=0.0, i0_percent=0.0)
case = "min" if np.isclose(c, 1.0) else "max"
fault_bus = bus1 if fault_node == 1 else bus2
pp.shortcircuit.calc_sc(net, case=case, ip=False, r_fault_ohm=0.0)
i_f_ka = net.res_bus_sc["ikss_ka"][fault_bus] # ka is actually amperes
return {"i_f": i_f_ka}
def calculate_vision(skss, c, fault_node, u_rated_transformer_node):
vision_results = {
(1000e6, 1.0, 1, 10e3): 5774,
(1000e6, 1.1, 1, 10e3): 5774,
(1000e6, 1.0, 2, 10e3): 52490,
(1000e6, 1.1, 2, 10e3): 52990,
(1000e6, 1.0, 2, 10.5e3): 55110,
(1000e6, 1.1, 2, 10.5e3): 55640,
}
return {"i_f": vision_results.get((skss, c, fault_node, u_rated_transformer_node), None)}
def calculate_hand(skss, c, fault_node, u_rated_transformer_node):
z_base_transformer_node = u_rated_transformer_node * u_rated_transformer_node / s_base
i_base_transformer_node = s_base / u_rated_transformer_node
z_k_pu = c / (skss / s_base)
if fault_node == 1:
i_f_pu = c * 1.0 / (z_k_pu * np.sqrt(3))
i_f = i_f_pu * i_base_source_node
return {
"i_f_iec": i_f,
}
elif fault_node == 2:
k = (u1 / u2) / (u_rated_source_node / u_rated_transformer_node)
z_k_lv_pu = z_k_pu * (1 / k / k)
z_k_lv_ohms = z_k_lv_pu * z_base_transformer_node
z_base_transformer_rating = u2 * u2 / sn
z_t_ohms = uk * z_base_transformer_rating
z_t_pu = z_t_ohms / z_base_transformer_node
z_total_pu = z_k_lv_pu + z_t_pu
z_total_ohms = z_k_lv_ohms + z_t_ohms
i_f_pu_pgm = (c * 1.0 / k) / (z_total_pu * np.sqrt(3))
i_f_pgm = i_f_pu_pgm * i_base_transformer_node
# Transformer impedance correction factor
if np.isclose(c, 1.1):
k_t_correction_factor = (0.95 * c) / ( 1 + 0.6 * uk)
z_total_ohms = z_k_lv_ohms + z_t_ohms * k_t_correction_factor
i_f_iec = c * u_rated_transformer_node/ (z_total_ohms * np.sqrt(3))
return {
"i_f_pgm": i_f_pgm,
"i_f_iec": i_f_iec,
}
def to_float(d):
return {k: (f"{float(v):.2f}" if isinstance(v, np.floating) else f"{v:.2f}") for k, v in d.items()}
def run():
cases = [
(1000e6, 1.0, 1, 10e3),
(1000e6, 1.1, 1, 10e3),
(1000e6, 1.0, 2, 10e3),
(1000e6, 1.0, 2, 10.5e3),
(1000e6, 1.1, 2, 10e3),
(1000e6, 1.1, 2, 10.5e3),
]
results = {}
for skss, c, fault_node, u_rated_transformer_node in cases:
results[(u_rated_transformer_node, c, fault_node)] = {
"pgm": to_float(calculate_pgm_engine(skss=skss, c=c, fault_node=fault_node, u_rated_transformer_node=u_rated_transformer_node)),
"hand": to_float(calculate_hand(skss=skss, c=c, fault_node=fault_node, u_rated_transformer_node=u_rated_transformer_node)),
"pandapower": to_float(calculate_pandapower(skss=skss, c=c, fault_node=fault_node, u_rated_transformer_node=u_rated_transformer_node)),
"vision": to_float(calculate_vision(skss=skss, c=c, fault_node=fault_node, u_rated_transformer_node=u_rated_transformer_node)),
}
print("\nCase header: (u_rated_transformer_node, c, fault_node)\n")
pprint(results, indent=4, sort_dicts=False, width=200)
if __name__ == "__main__":
run()gives output
Case header: (u_rated_transformer_node, c, fault_node)
{ (10000.0, 1.0, 1): { 'pgm': {'i_f': '5773.50'},
'hand': {'i_f_iec': '5773.50'},
'pandapower': {'i_f': '5773.50'},
'vision': {'i_f': 5774}},
(10000.0, 1.1, 1): { 'pgm': {'i_f': '5773.50'},
'hand': {'i_f_iec': '5773.50'},
'pandapower': {'i_f': '5773.50'},
'vision': {'i_f': 5774}},
(10000.0, 1.0, 2): { 'pgm': {'i_f': '52486.39'},
'hand': {'i_f_pgm': '52486.39', 'i_f_iec': '52486.39'},
'pandapower': {'i_f': '52554.00'},
'vision': {'i_f': 52490}},
(10500.0, 1.0, 2): { 'pgm': {'i_f': '52486.39'},
'hand': {'i_f_pgm': '52486.39', 'i_f_iec': '55110.71'},
'pandapower': {'i_f': '55181.70'},
'vision': {'i_f': 55110}},
(10000.0, 1.1, 2): { 'pgm': {'i_f': '52923.77'},
'hand': {'i_f_pgm': '52923.77', 'i_f_iec': '52986.26'},
'pandapower': {'i_f': '52986.26'},
'vision': {'i_f': 52990}},
(10500.0, 1.1, 2): { 'pgm': {'i_f': '52923.77'},
'hand': {'i_f_pgm': '52923.77', 'i_f_iec': '55635.57'},
'pandapower': {'i_f': '55635.57'},
'vision': {'i_f': 55640}}}
Note:
- Pandapower minimum results include the transformer impedance correction as well (I see conflicting information regarding this in the internet, need to verify which is correct). If we compare it with a hypothetical hand calculation including c=1.1 impedance correction factor, we get the exact same result. Maybe this is a bug at their end, I shall raise an issue there after we resolve our discussion here.
- Vision results are rounded off.
IEC standard specifies usage of rated voltage of the node to do short circuit calculations. This is because it makes no assumption on actual voltage but assumes that nominally it would be 1.0 p.u. for a bus and at extreme, 1.1.
PGM is not 100% conformant with IEC 60909.
We do not calculate short circuit currents using rated voltage of fault nodes:
c * u_rated. But we usec * u_rated * k_eqwhich includes any transformation ratiok_eqinvolved from source to the fault node. Thiskfactor is statically present in the grid and not dependent on load conditions hence it is present in exact form in the short circuit situation too.The design decision was that since our approach would be more accurate compared to the standard, we do not have to stick with certain assumptions.
The other decisions like inclusion of magnetisation current, tap positions and capcitance in the short circuit calculations are easily managed on the user side by setting i0, c1, c0, tap_pos, tap_nom as 0. In the case of voltage, determining
k_eqespecially with multiple sources becomes more complex to bypass.Many of MV grids have varied transformer secondary ratings compared to a fixed 10kV node voltage rating, hence this issue concerns those cases.
As expected, we also deviate from other calculation engines in this regard (See in the script below). Hence our validation does not become straightforward.
Given these points, let us revisit this decision again.
Comparison script
Reproduces PGM calculations, calculations inside PGM replicated by hand, hand calculations via IEC standard, vision (Calculated externally) and pandapower. (TLDR to relevant functions or output directly.)
The parameters are unrealistic but suitable for understanding the calculation.
gives output
Note: