From c8a37e8a75f449a2095abbe93c0b085ae9a7a4ac Mon Sep 17 00:00:00 2001 From: Mike Vogt Date: Mon, 31 Mar 2025 10:21:06 +0200 Subject: [PATCH 1/2] first push of revised harmonic calculation code. --- pandapower/harmonics/__init__.py | 0 pandapower/harmonics/balanced.py | 151 +++++++++ .../harmonics/harmonic_impedance_creator.py | 194 +++++++++++ pandapower/harmonics/unbalanced.py | 302 ++++++++++++++++++ pandapower/test/harmonics/Test_Case_1.py | 47 +++ pandapower/test/harmonics/Test_Case_2.py | 39 +++ pandapower/test/harmonics/Test_Case_3.py | 46 +++ pandapower/test/harmonics/Test_Case_4.py | 87 +++++ pandapower/test/harmonics/Test_Case_5.py | 78 +++++ pandapower/test/harmonics/Test_Case_6.py | 82 +++++ pandapower/test/harmonics/__init__.py | 0 .../test/test_files/harmonics/CIGRE_LV.xlsx | Bin 0 -> 12363 bytes .../test_files/harmonics/Test_1f_Loads.xlsx | Bin 0 -> 8629 bytes .../test_files/harmonics/Test_1f_lines.xlsx | Bin 0 -> 9952 bytes 14 files changed, 1026 insertions(+) create mode 100644 pandapower/harmonics/__init__.py create mode 100644 pandapower/harmonics/balanced.py create mode 100644 pandapower/harmonics/harmonic_impedance_creator.py create mode 100644 pandapower/harmonics/unbalanced.py create mode 100644 pandapower/test/harmonics/Test_Case_1.py create mode 100644 pandapower/test/harmonics/Test_Case_2.py create mode 100644 pandapower/test/harmonics/Test_Case_3.py create mode 100644 pandapower/test/harmonics/Test_Case_4.py create mode 100644 pandapower/test/harmonics/Test_Case_5.py create mode 100644 pandapower/test/harmonics/Test_Case_6.py create mode 100644 pandapower/test/harmonics/__init__.py create mode 100644 pandapower/test/test_files/harmonics/CIGRE_LV.xlsx create mode 100644 pandapower/test/test_files/harmonics/Test_1f_Loads.xlsx create mode 100644 pandapower/test/test_files/harmonics/Test_1f_lines.xlsx diff --git a/pandapower/harmonics/__init__.py b/pandapower/harmonics/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pandapower/harmonics/balanced.py b/pandapower/harmonics/balanced.py new file mode 100644 index 0000000000..d4544f6ae5 --- /dev/null +++ b/pandapower/harmonics/balanced.py @@ -0,0 +1,151 @@ +import math +import cmath +import numpy as np +import pandapower as pp +import pandapower.harmonics.harmonic_impedance_creator as hic + + +# formatting harmonic voltages in table and calculation of THD +def balanced_thd_voltage(net, harmonics, har, har_angle, analysis_type): + harmonics_voltage_0, harmonics_voltage = balanced_harmonic_current_voltage(net, harmonics, har, + har_angle, analysis_type) + + # Nodes need to be sorted 0, 1, 2, 3, 4... Node with index 0 needs to be referent node (External grid is connected to this node) + + thd = [] + for i in range(0, np.shape(harmonics_voltage)[0]): + sum_thd = 0 + for j in range(0, np.shape(harmonics_voltage)[1]): + sum_thd += harmonics_voltage[i, j] ** 2 + thd.append(sum_thd) + + for i in range(0, len(thd)): + thd[i] = math.sqrt(thd[i]) / (net.res_bus.vm_pu[i + 1]) * 100 + + sum_thd_0 = 0 + + for i in range(0, len(harmonics_voltage_0)): + sum_thd_0 += harmonics_voltage_0[i] ** 2 + + thd.insert(0, math.sqrt(sum_thd_0) / net.res_bus.vm_pu[0] * 100) + + harmonics_voltage_res = np.zeros([int(len(net.bus.name)), len(har)], dtype=float) + + for i in range(0, np.shape(harmonics_voltage)[0]): + for j in range(0, np.shape(harmonics_voltage)[1]): + harmonics_voltage_res[i + 1, j] = harmonics_voltage[i, j] * 100 + + for i in range(0, len(harmonics_voltage_0)): + harmonics_voltage_res[0, i] = harmonics_voltage_0[i] * 100 + + return thd, harmonics_voltage_res + + +# calculation of harmonic voltages from harmonic currents +# at the moment it is possible to define only one harmonic patter which is same for every harmonic source +def balanced_harmonic_current_voltage(net, harmonics, har, har_angle, analysis_type): + delta_harmonics_voltage = np.zeros([len(net.bus.name) - 1, len(har)], dtype=complex) + harmonics_voltage = np.zeros([len(net.bus.name) - 1, len(har)], dtype=float) + harmonic_cur_val = [] + harmonic_cur_ang = [] + u_harmonics_0 = [] + u_harmonics_0_ang = [] + harmonic_cur_0 = [] + harmonic_cur_0_ang = [] + + har_matrices, har_ext_matrix = hic.harmonic_imp_creator(net, harmonics, analysis_type) + + for h in range(0, len(harmonics)): + mat_z = har_matrices[h] + z_ext = har_ext_matrix[h] + + if harmonics[h] == 1: + pp.runpp(net) + + current = [] + current_0 = [] + + s = complex(net.res_bus.p_mw[net.ext_grid.bus[0]], net.res_bus.q_mvar[net.ext_grid.bus[0]]) + + current_0.append(np.conjugate(s / (math.sqrt(3) * (cmath.rect(net.res_bus.vm_pu[net.ext_grid.bus[0]], \ + net.res_bus.va_degree[ + net.ext_grid.bus[0]]))))) + + for i in net.res_bus.index: + connected = 0 + + if i != net.ext_grid.bus[0]: + for j in range(0, len(net.load.index)): + + if net.load.bus[j] == net.res_bus.index[i]: + connected = 1 + + s = complex(net.res_bus.p_mw[i], net.res_bus.q_mvar[i]) + + if connected == 1: + current.append(np.conjugate(s \ + / (math.sqrt(3) * ( + cmath.rect(net.res_bus.vm_pu[i], net.res_bus.va_degree[i]))))) + else: + current.append(0 + 0j) + else: + + harmonics_current = [] + harmonics_angle = [] + + for i in range(0, len(net.bus.name)): + harmonics_current.append(har[h - 1]) + # Only positive sequence system is considered. harmonics[h]*har_a_angle[h-1] + 240 can be appended + # but it does not change the magnitude of the voltage, only the angle. + if analysis_type == 'balanced_positive': + + if harmonics[h] % 3 == 0: + harmonics_angle.append(harmonics[h] * har_angle[h - 1] + 240) + elif harmonics[h] % 3 == 1: + harmonics_angle.append(harmonics[h] * har_angle[h - 1] + 240) + elif harmonics[h] % 3 == 2: + harmonics_angle.append(harmonics[h] * har_angle[h - 1] + 240) + + elif analysis_type == 'balanced_all': + + if harmonics[h] % 3 == 0: + harmonics_angle.append(harmonics[h] * har_angle[h - 1]) + elif harmonics[h] % 3 == 1: + harmonics_angle.append(harmonics[h] * har_angle[h - 1] + 240) + elif harmonics[h] % 3 == 2: + harmonics_angle.append(harmonics[h] * har_angle[h - 1] + 120) + + har_cur = [] + + for i in range(0, len(current)): + har_cur.append( + -abs(current[i]) * cmath.rect(harmonics_current[i] / 100, harmonics_angle[i] * cmath.pi / 180)) + harmonic_cur_val.append(abs(har_cur[i])) + harmonic_cur_ang.append(cmath.phase(har_cur[i]) * 180 / cmath.pi) + + sum_a = 0 + + for a in range(0, len(har_cur)): + sum_a += har_cur[a] + + har_cur_0 = [] + har_cur_0.append(sum_a) + + harmonic_cur_0.append(abs(sum_a)) + harmonic_cur_0_ang.append(cmath.phase(sum_a) * 180 / cmath.pi) + + u_har_0 = (z_ext * har_cur_0[0] * math.sqrt(3)) + + u_harmonics_0.append(abs(u_har_0)) + u_harmonics_0_ang.append(cmath.phase(u_har_0) * 180 / cmath.pi) + + delta_har_vol = np.matmul(mat_z, har_cur) * math.sqrt(3) + + for i in range(0, np.shape(delta_harmonics_voltage)[0]): + delta_harmonics_voltage[i, h - 1] = (np.transpose(delta_har_vol)[i]) + + for i in range(0, np.shape(harmonics_voltage)[0]): + harmonics_voltage[i, h - 1] = abs(u_har_0 + delta_harmonics_voltage[i, h - 1]) + + return u_harmonics_0, harmonics_voltage + diff --git a/pandapower/harmonics/harmonic_impedance_creator.py b/pandapower/harmonics/harmonic_impedance_creator.py new file mode 100644 index 0000000000..858c9e99ff --- /dev/null +++ b/pandapower/harmonics/harmonic_impedance_creator.py @@ -0,0 +1,194 @@ +#Script for creating impedance matrix for higher order frequencies + +import pandapower as pp +import numpy as np +import math +import cmath +from pandapower.pd2ppc import _pd2ppc +from pandapower.pd2ppc_zero import _pd2ppc_zero + +def harmonic_imp_creator(net, harmon_order, analysis_type): + harmonics = harmon_order + + a1 = cmath.rect(1, 2/3*cmath.pi) + a2 = cmath.rect(1, 4/3*cmath.pi) + + matrix_A = np.matrix([[1, 1, 1], [1, a2, a1], [1, a1, a2]]) + + harmonic_matrices = [] + harmonic_ext_matrices = [] + + #In case of unbalanced harmonic analyses, values for zero sequence system need to be defined. + #If they are unknown, they are assumed to be same as values for positive sequence system. + if 'r0_ohm_per_km' not in net.line: + + net.line['r0_ohm_per_km'] = net.line['r_ohm_per_km'] + + if 'x0_ohm_per_km' not in net.line: + net.line['x0_ohm_per_km'] = net.line['x_ohm_per_km'] + + if 'c0_nf_per_km' not in net.line: + net.line['c0_nf_per_km'] = net.line['c_nf_per_km'] + + #Defining parameters needed for LF calculations + + net["_options"] = {"mode":"pf", "check_connectivity":True, "calculate_voltage_angles":True, "init_vm_pu":True, + "init_va_degree":True, "consider_line_temperature":False, "voltage_depend_loads":False, + "trafo3w_losses":"hv", "neglect_open_switch_branches":False, "p_lim_default":1000000, + "delta":0, "q_lim_default":1000000, "trafo_model":"t", "distributed_slack": False, "tdpf": False} + net["_isolated_buses"] = [] + net["_is_elements"] = None + + u_base = net.bus.vn_kv[net.ext_grid.bus[0]] + + #Calculation of external grid's impedance and the network's impedance matrix for every harmonic order + #Already developed pandapower functionalities are used for matrices creation + for h in range(0, len(harmonics)): + + ppc_0, ppci_0 = _pd2ppc_zero(net, None) + + ybus_0, yf_0, yt_0 = pp.pypower.makeYbus.makeYbus(ppci_0["baseMVA"], ppci_0["bus"], ppci_0["branch"]) + + ppc_1, ppci_1 = _pd2ppc(net, 1) + + ybus_1, yf_1, yt_1 = pp.pypower.makeYbus.makeYbus(ppci_1["baseMVA"], ppci_1["bus"], ppci_1["branch"]) + + ppc_2, ppci_2 = _pd2ppc(net, 2) + + ybus_2, yf_2, yt_2 = pp.pypower.makeYbus.makeYbus(ppci_2["baseMVA"], ppci_2["bus"], ppci_2["branch"]) + + #Full sequence admittance matrices + zero_full = ybus_0.todense() + pos_full = ybus_1.todense() + neg_full = ybus_2.todense() + + #Removing referent node related row and column of each matrix + #It is assumed that first row and column are referent node related + zero = np.delete(np.delete(zero_full, 0, 0), 0, 1) + pos = np.delete(np.delete(pos_full, 0, 0), 0, 1) + neg = np.delete(np.delete(neg_full, 0, 0), 0, 1) + + #External grid impedance calculation + + z_pos_ohm = u_base**2/net.ext_grid.s_sc_max_mva + delta_pos = math.atan(1/net.ext_grid.rx_max) + r_pos = z_pos_ohm*math.cos(delta_pos) + x_pos = z_pos_ohm*math.sin(delta_pos)*harmonics[h] + z_pos = complex(r_pos, x_pos) + + z_zer_ohm = z_pos_ohm*net.ext_grid.x0x_max + delta_zer = math.atan(1/net.ext_grid.r0x0_max) + r_zer = z_zer_ohm*math.cos(delta_zer) + x_zer = z_zer_ohm*math.sin(delta_zer)*harmonics[h] + z_zer = complex(r_zer, x_zer) + + #Depenent on the analysis type, impedances are created + if analysis_type == 'unbalanced': + + phase_mat_y = np.zeros([3*(np.shape(zero)[0]), 3*(np.shape(zero)[1])], dtype = complex) + + z_ext_0 = np.zeros([3,3], dtype = complex) + + #Z012 matrix is created for the referent node + z_ext_0[0, 0] = z_zer/(u_base**2/3) + z_ext_0[1, 1] = z_pos/(u_base**2/3) + z_ext_0[2, 2] = z_pos/(u_base**2/3) + + #Transformation from the sequence to the phase system + z_ext_abc = np.matmul(np.matmul(np.linalg.inv(matrix_A), z_ext_0), matrix_A) + + for i in range(0, np.shape(zero)[0]): + for j in range(0, np.shape(zero)[1]): + #Y012 is created for every other element of Y0,1,2 matrices + y_012 = np.zeros([3,3], dtype = complex) + + z = zero[i, j] + p = pos[i,j] + n = neg[i,j] + + y_012[0,0] = z + y_012[1,1] = p + y_012[2,2] = n + + #Transformation from the sequence to the phase system + y_abc = np.matmul(np.matmul(np.linalg.inv(matrix_A), y_012), matrix_A) + + for row in range (0, 3): + for col in range(0, 3): + phase_mat_y[3*i+row, 3*j+col] = y_abc[row, col] + #Z matrix as an inverse of the Y matrix + phase_mat_z = np.linalg.inv(phase_mat_y)*3 + + #Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X + for r in range(0, np.shape(phase_mat_z)[0]): + for c in range(0, np.shape(phase_mat_z)[1]): + aux = phase_mat_z[r,c] + res = np.real(aux) + reac = np.imag(aux) * harmonics[h] + imp = complex(res, reac) + phase_mat_z[r,c] = imp + + #Appending the referent node and network impedances + harmonic_matrices.append(phase_mat_z) + harmonic_ext_matrices.append(z_ext_abc) + + elif analysis_type == 'balanced_positive': + + #Only positive sequence system is observed + z_pos_net = np.linalg.inv(pos) + + #Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X + for r in range(0, np.shape(z_pos_net)[0]): + for c in range(0, np.shape(z_pos_net)[1]): + aux = z_pos_net[r,c] + res = np.real(aux) + reac = np.imag(aux) * harmonics[h] + imp = complex(res, reac) + z_pos_net[r,c] = imp + + #Appending the referent node and network impedances + harmonic_matrices.append(z_pos_net) + harmonic_ext_matrices.append(z_pos/(u_base**2)) + + elif analysis_type == 'balanced_all': + + #Balanced network but all sequence systems are used + #Dependent on the harmonic order, zero, positive or negative sequence impedance are considered + z_pos_net = np.linalg.inv(pos) + z_neg_net = np.linalg.inv(neg) + z_zero_net = np.linalg.inv(zero) + + #Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X + for r in range(0, np.shape(z_pos_net)[0]): + for c in range(0, np.shape(z_pos_net)[1]): + aux_p = z_pos_net[r,c] + res_p = np.real(aux_p) + reac_p = np.imag(aux_p) * harmonics[h] + imp_p = complex(res_p, reac_p) + z_pos_net[r,c] = imp_p + + aux_n = z_neg_net[r,c] + res_n = np.real(aux_n) + reac_n = np.imag(aux_n) * harmonics[h] + imp_n = complex(res_n, reac_n) + z_neg_net[r,c] = imp_n + + aux_z = z_zero_net[r,c] + res_z = np.real(aux_z) + reac_z = np.imag(aux_z) * harmonics[h] + imp_z = complex(res_z, reac_z) + z_zero_net[r,c] = imp_z + + #Appending the referent node and network impedances + if harmonics[h] % 3 == 0: + harmonic_matrices.append(z_zero_net) + harmonic_ext_matrices.append(z_zer/(u_base**2)) + elif harmonics[h] % 3 == 1: + harmonic_matrices.append(z_pos_net) + harmonic_ext_matrices.append(z_pos/(u_base**2)) + elif harmonics[h] % 3 == 2: + harmonic_matrices.append(z_neg_net) + harmonic_ext_matrices.append(z_pos/(u_base**2)) + + return(harmonic_matrices, harmonic_ext_matrices) + \ No newline at end of file diff --git a/pandapower/harmonics/unbalanced.py b/pandapower/harmonics/unbalanced.py new file mode 100644 index 0000000000..c2a99ca423 --- /dev/null +++ b/pandapower/harmonics/unbalanced.py @@ -0,0 +1,302 @@ +import math +import cmath +import numpy as np +import pandapower as pp +import pandapower.harmonics.harmonic_impedance_creator as hic +#import pandapower.harmonics.unbalanced.voltages_currents_calculation as har_vol_calc + + +# formatting harmonic voltages in table and calculation of THD +# Nodes need to be sorted 0, 1, 2, 3, 4... Node with index 0 needs to be referent node (External grid is connected to this node) + +def unbalanced_thd_voltage(network, harmonics, har_a, har_a_angle, har_b, har_b_angle, har_c, har_c_angle, + har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle, + analysis_type): + harmonics_voltage_0, harmonics_voltage, harmonics_current_0, harmonics_current = \ + unbalanced_harmonic_current_voltage(network, harmonics, + har_a, har_a_angle, har_b, har_b_angle, har_c, har_c_angle, + har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, + har_c_lc_angle, + analysis_type) + + thd_a = [] + thd_b = [] + thd_c = [] + + for i in range(0, np.shape(harmonics_voltage)[0], 3): + sum_thd = 0 + for j in range(0, np.shape(harmonics_voltage)[1]): + sum_thd += harmonics_voltage[i, j] ** 2 + thd_a.append(sum_thd) + + for i in range(1, np.shape(harmonics_voltage)[0], 3): + sum_thd = 0 + for j in range(0, np.shape(harmonics_voltage)[1]): + sum_thd += harmonics_voltage[i, j] ** 2 + thd_b.append(sum_thd) + + for i in range(2, np.shape(harmonics_voltage)[0], 3): + sum_thd = 0 + for j in range(0, np.shape(harmonics_voltage)[1]): + sum_thd += harmonics_voltage[i, j] ** 2 + thd_c.append(sum_thd) + + for i in range(0, len(thd_a)): + thd_a[i] = math.sqrt(thd_a[i]) / (network.res_bus_3ph.vm_a_pu[i + 1]) * 100 + thd_b[i] = math.sqrt(thd_b[i]) / (network.res_bus_3ph.vm_b_pu[i + 1]) * 100 + thd_c[i] = math.sqrt(thd_c[i]) / (network.res_bus_3ph.vm_c_pu[i + 1]) * 100 + + sum_thd_a_0 = 0 + sum_thd_b_0 = 0 + sum_thd_c_0 = 0 + + for i in range(0, np.shape(harmonics_voltage_0)[1]): + sum_thd_a_0 += harmonics_voltage_0[0, i] ** 2 + sum_thd_b_0 += harmonics_voltage_0[1, i] ** 2 + sum_thd_c_0 += harmonics_voltage_0[2, i] ** 2 + + thd_a.insert(0, math.sqrt(sum_thd_a_0) / network.res_bus_3ph.vm_a_pu[0] * 100) + thd_b.insert(0, math.sqrt(sum_thd_b_0) / network.res_bus_3ph.vm_b_pu[0] * 100) + thd_c.insert(0, math.sqrt(sum_thd_c_0) / network.res_bus_3ph.vm_c_pu[0] * 100) + + harmonics_voltage_a = np.zeros([int(len(network.bus.name) * 3 / 3), len(har_a)], dtype=float) + harmonics_voltage_b = np.zeros([int(len(network.bus.name) * 3 / 3), len(har_a)], dtype=float) + harmonics_voltage_c = np.zeros([int(len(network.bus.name) * 3 / 3), len(har_a)], dtype=float) + # Harmonic voltage (percentage) of each phase and each harmoni + for i in range(0, np.shape(harmonics_voltage_a)[0] - 1): + for j in range(0, np.shape(harmonics_voltage)[1]): + harmonics_voltage_a[i + 1, j] = harmonics_voltage[3 * i, j] * 100 * (1 / network.res_bus_3ph.vm_a_pu[i + 1]) + harmonics_voltage_b[i + 1, j] = harmonics_voltage[3 * i + 1, j] * 100 * ( + 1 / network.res_bus_3ph.vm_b_pu[i + 1]) + harmonics_voltage_c[i + 1, j] = harmonics_voltage[3 * i + 2, j] * 100 * ( + 1 / network.res_bus_3ph.vm_c_pu[i + 1]) + + for i in range(0, np.shape(harmonics_voltage_0)[1]): + harmonics_voltage_a[0, i] = harmonics_voltage_0[0, i] * 100 * (1 / network.res_bus_3ph.vm_a_pu[0]) + harmonics_voltage_b[0, i] = harmonics_voltage_0[1, i] * 100 * (1 / network.res_bus_3ph.vm_b_pu[0]) + harmonics_voltage_c[0, i] = harmonics_voltage_0[2, i] * 100 * (1 / network.res_bus_3ph.vm_c_pu[0]) + + return thd_a, thd_b, thd_c, harmonics_voltage_a, harmonics_voltage_b, harmonics_voltage_c + + +def unbalanced_harmonic_current_voltage(network, harmonics, har_a, har_a_angle, har_b, har_b_angle, har_c, har_c_angle, \ + har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle, \ + analysis_type): + delta_harmonics_voltage = np.zeros([len(network.bus.name) * 3 - 3, len(har_a)], dtype=complex) + harmonics_voltage = np.zeros([len(network.bus.name) * 3 - 3, len(har_a)], dtype=float) + harmonics_current = np.zeros([len(network.bus.name) * 3 - 3, len(har_a)], dtype=float) + u_harmonics_0 = np.zeros([3, len(har_a)], dtype=complex) + i_harmonics_0 = np.zeros([3, len(har_a)], dtype=float) + + har_matrices, har_ext_matrix = hic.harmonic_imp_creator(network, harmonics, analysis_type) + + s_base = network.sn_mva * 1e6 + + for h in range(0, len(harmonics)): + phase_mat_z = har_matrices[h] + z_ext_abc = har_ext_matrix[h] + + if harmonics[h] == 1: + + pp.runpp_3ph(network) + + current = [] + current_res = [] + current_pv = [] + current_0 = [] + + # p.u. values + s_a = complex(network.res_bus_3ph.p_a_mw[network.ext_grid.bus[0]] * 1000000 / s_base, + network.res_bus_3ph.q_a_mvar[network.ext_grid.bus[0]] * 1000000 / s_base) + s_b = complex(network.res_bus_3ph.p_b_mw[network.ext_grid.bus[0]] * 1000000 / s_base, + network.res_bus_3ph.q_b_mvar[network.ext_grid.bus[0]] * 1000000 / s_base) + s_c = complex(network.res_bus_3ph.p_c_mw[network.ext_grid.bus[0]] * 1000000 / s_base, + network.res_bus_3ph.q_c_mvar[network.ext_grid.bus[0]] * 1000000 / s_base) + + current_0.append(np.conjugate(s_a / (cmath.rect(network.res_bus_3ph.vm_a_pu[network.ext_grid.bus[0]], + network.res_bus_3ph.va_a_degree[ + network.ext_grid.bus[0]] * cmath.pi / 180)))) + current_0.append(np.conjugate(s_b / (cmath.rect(network.res_bus_3ph.vm_b_pu[network.ext_grid.bus[0]], + network.res_bus_3ph.va_b_degree[ + network.ext_grid.bus[0]] * cmath.pi / 180)))) + current_0.append(np.conjugate(s_c / (cmath.rect(network.res_bus_3ph.vm_c_pu[network.ext_grid.bus[0]], + network.res_bus_3ph.va_c_degree[ + network.ext_grid.bus[0]] * cmath.pi / 180)))) + + for i in network.res_bus_3ph.index: + connected = 0 + + if i != network.ext_grid.bus[0]: + for j in range(0, len(network.asymmetric_load.index), 2): + # The assumption is that two harmonic sources are connected to the node + # It can be meodified, additional sources can be added, or number of sources at every node can be reduced + if network.asymmetric_load.bus[j] == i: + connected = 1 # an asymmetric load is connected to the observed node + + ##Residential + s_a_res = complex(network.asymmetric_load.p_a_mw[j] * 1000000 / s_base, + network.asymmetric_load.q_a_mvar[j] * 1000000 / s_base) + s_b_res = complex(network.asymmetric_load.p_b_mw[j] * 1000000 / s_base, + network.asymmetric_load.q_b_mvar[j] * 1000000 / s_base) + s_c_res = complex(network.asymmetric_load.p_c_mw[j] * 1000000 / s_base, + network.asymmetric_load.q_c_mvar[j] * 1000000 / s_base) + + # PV + s_a_pv = complex(network.asymmetric_load.p_a_mw[j + 1] * 1000000 / s_base, + network.asymmetric_load.q_a_mvar[j + 1] * 1000000 / s_base) + s_b_pv = complex(network.asymmetric_load.p_b_mw[j + 1] * 1000000 / s_base, + network.asymmetric_load.q_b_mvar[j + 1] * 1000000 / s_base) + s_c_pv = complex(network.asymmetric_load.p_c_mw[j + 1] * 1000000 / s_base, + network.asymmetric_load.q_c_mvar[j + 1] * 1000000 / s_base) + + break + + if connected == 1: + # If the node is a node where a load/generator is connected we calculate the currents + # Else currents at that node are 0 + current_res.append(np.conjugate(s_a_res / (cmath.rect(network.res_bus_3ph.vm_a_pu[i], + network.res_bus_3ph.va_a_degree[ + i] * cmath.pi / 180)))) + current_res.append(np.conjugate(s_b_res / (cmath.rect(network.res_bus_3ph.vm_b_pu[i], + network.res_bus_3ph.va_b_degree[ + i] * cmath.pi / 180)))) + current_res.append(np.conjugate(s_c_res / (cmath.rect(network.res_bus_3ph.vm_c_pu[i], + network.res_bus_3ph.va_c_degree[ + i] * cmath.pi / 180)))) + + current_pv.append(np.conjugate(s_a_pv / (cmath.rect(network.res_bus_3ph.vm_a_pu[i], + network.res_bus_3ph.va_a_degree[ + i] * cmath.pi / 180)))) + current_pv.append(np.conjugate(s_b_pv / (cmath.rect(network.res_bus_3ph.vm_b_pu[i], + network.res_bus_3ph.va_b_degree[ + i] * cmath.pi / 180)))) + current_pv.append(np.conjugate(s_c_pv / (cmath.rect(network.res_bus_3ph.vm_c_pu[i], + network.res_bus_3ph.va_c_degree[ + i] * cmath.pi / 180)))) + + else: + current_res.append(0 + 0j) + current_res.append(0 + 0j) + current_res.append(0 + 0j) + + current_pv.append(0 + 0j) + current_pv.append(0 + 0j) + current_pv.append(0 + 0j) + + # p.u. values of current in each node. Must be equal to s_res + s_pv + s_a = complex(network.res_bus_3ph.p_a_mw[i] * 1000000 / s_base, + network.res_bus_3ph.q_a_mvar[i] * 1000000 / s_base) + s_b = complex(network.res_bus_3ph.p_b_mw[i] * 1000000 / s_base, + network.res_bus_3ph.q_b_mvar[i] * 1000000 / s_base) + s_c = complex(network.res_bus_3ph.p_c_mw[i] * 1000000 / s_base, + network.res_bus_3ph.q_c_mvar[i] * 1000000 / s_base) + + current.append(np.conjugate(s_a / (cmath.rect(network.res_bus_3ph.vm_a_pu[i], + network.res_bus_3ph.va_a_degree[ + i] * cmath.pi / 180)))) + current.append(np.conjugate(s_b / (cmath.rect(network.res_bus_3ph.vm_b_pu[i], + network.res_bus_3ph.va_b_degree[ + i] * cmath.pi / 180)))) + current.append(np.conjugate(s_c / (cmath.rect(network.res_bus_3ph.vm_c_pu[i], + network.res_bus_3ph.va_c_degree[ + i] * cmath.pi / 180)))) + + else: + harmonics_current_res = [] + harmonics_angle_res = [] + + harmonics_current_pv = [] + harmonics_angle_pv = [] + + for i in range(0, len(network.bus.name)): + # Residential household devices harmonics + harmonics_current_res.append(har_a[h - 1]) + harmonics_current_res.append(har_b[h - 1]) + harmonics_current_res.append(har_c[h - 1]) + + # PV harmonics + harmonics_current_pv.append(har_a_lc[h - 1]) + harmonics_current_pv.append(har_b_lc[h - 1]) + harmonics_current_pv.append(har_c_lc[h - 1]) + + # Harmonic angles are defined. The harmonic angles are not same for 3rd, 5th, 7th etc. + + if harmonics[h] % 3 == 0: + harmonics_angle_res.append(harmonics[h] * har_a_angle[h - 1]) + harmonics_angle_res.append(harmonics[h] * har_b_angle[h - 1]) + harmonics_angle_res.append(harmonics[h] * har_c_angle[h - 1]) + + harmonics_angle_pv.append(harmonics[h] * har_a_lc_angle[h - 1]) + harmonics_angle_pv.append(harmonics[h] * har_b_lc_angle[h - 1]) + harmonics_angle_pv.append(harmonics[h] * har_c_lc_angle[h - 1]) + elif harmonics[h] % 3 == 1: + harmonics_angle_res.append(harmonics[h] * har_a_angle[h - 1]) + harmonics_angle_res.append(harmonics[h] * har_b_angle[h - 1] + 240) + harmonics_angle_res.append(harmonics[h] * har_c_angle[h - 1] + 120) + + harmonics_angle_pv.append(harmonics[h] * har_a_lc_angle[h - 1]) + harmonics_angle_pv.append(harmonics[h] * har_b_lc_angle[h - 1] + 240) + harmonics_angle_pv.append(harmonics[h] * har_c_lc_angle[h - 1] + 120) + elif harmonics[h] % 3 == 2: + harmonics_angle_res.append(harmonics[h] * har_a_angle[h - 1]) + harmonics_angle_res.append(harmonics[h] * har_b_angle[h - 1] + 120) + harmonics_angle_res.append(harmonics[h] * har_c_angle[h - 1] + 240) + + harmonics_angle_pv.append(harmonics[h] * har_a_lc_angle[h - 1]) + harmonics_angle_pv.append(harmonics[h] * har_b_lc_angle[h - 1] + 120) + harmonics_angle_pv.append(harmonics[h] * har_c_lc_angle[h - 1] + 240) + + har_cur_pv = [] + har_cur_res = [] + har_cur = [] + + for i in range(0, len(current)): + # Only the absolute value of the fundamnetal harmonic is taken + # since in the most researches and measurments the angle of the fundamental harmonic is taken to be zero + # We caclulate harmonic currents of hosuehold devices, PVs, and we create harmonic current at the node + # as the sum of currents of all harmonic sources at the node + har_cur_res.append(-abs(current_res[i]) * cmath.rect(harmonics_current_res[i] / 100, + harmonics_angle_res[i] * cmath.pi / 180)) + har_cur_pv.append(-abs(current_pv[i]) * cmath.rect(harmonics_current_pv[i] / 100, + harmonics_angle_pv[i] * cmath.pi / 180)) + har_cur.append(har_cur_pv[i] + har_cur_res[i]) + + sum_a = 0 + sum_b = 0 + sum_c = 0 + + # At the first node, current is equal to the sum of harmonic currents at all other nodes + for a in range(0, len(har_cur), 3): + sum_a += har_cur[a] + sum_b += har_cur[a + 1] + sum_c += har_cur[a + 2] + + har_cur_0 = [] + + har_cur_0.append(sum_a) + har_cur_0.append(sum_b) + har_cur_0.append(sum_c) + + # Calculation of the harmonic current of the first node + u_har_0 = np.matmul(z_ext_abc, har_cur_0) + + for i in range(0, np.shape(u_har_0)[1]): + u_harmonics_0[i, h - 1] = u_har_0[0, i] + i_harmonics_0[i, h - 1] = abs(har_cur_0[i]) + + # Calculating the harmonics voltage drop, the difference between the harmonic of slack and all other nodes + delta_har_vol = np.matmul(phase_mat_z, har_cur) + + for i in range(0, np.shape(delta_harmonics_voltage)[0]): + delta_harmonics_voltage[i, h - 1] = delta_har_vol[i] + harmonics_current[i, h - 1] = abs(har_cur[i]) + + # Calculating harmonic voltage from the harmonic of a slack node and differences + for i in range(0, np.shape(delta_harmonics_voltage)[0], 3): + harmonics_voltage[i, h - 1] = abs(u_harmonics_0[0, h - 1] + delta_harmonics_voltage[i, h - 1]) + harmonics_voltage[i + 1, h - 1] = abs(u_harmonics_0[1, h - 1] + delta_harmonics_voltage[i + 1, h - 1]) + harmonics_voltage[i + 2, h - 1] = abs(u_harmonics_0[2, h - 1] + delta_harmonics_voltage[i + 2, h - 1]) + + harmonics_voltage_0 = abs(u_harmonics_0) + + return harmonics_voltage_0, harmonics_voltage, i_harmonics_0, harmonics_current \ No newline at end of file diff --git a/pandapower/test/harmonics/Test_Case_1.py b/pandapower/test/harmonics/Test_Case_1.py new file mode 100644 index 0000000000..efd9e5951a --- /dev/null +++ b/pandapower/test/harmonics/Test_Case_1.py @@ -0,0 +1,47 @@ +import os +import pytest +import pandas as pd +import pandapower as pp +from pandapower.harmonics.balanced import balanced_thd_voltage +from pandapower.plotting.simple_plot import simple_plot + +def test_harmonics_case_1(): + path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") + + nodes = pd.read_excel(os.path.join(path, 'Test_1f_Loads.xlsx')) + lines = pd.read_excel(os.path.join(path,'Test_1f_lines.xlsx')) + + har_a = [11,10,9,8,7,6,5] + har_a_angle = [250,240,230,220,210,200,190] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + # harmonics = [1, 3] + + net = pp.create_empty_network() + + for i in range(0, len(nodes['Naziv'])): + pp.create_bus(net, 10, nodes['Naziv'][i], int(nodes['ID'][i])) + + pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 10, s_sc_min_mva = 10, rx_max=5, rx_min=0, r0x0_max=5, x0x_max=3) + + for i in range (0, len(lines['Length'])): + pp.create_line_from_parameters(net, + lines['Start Node'][i], + lines['End Node'][i], + lines['Length'][i], + lines['R1[ohm/km]'][i], + lines['X1[ohm/km]'][i], + 1e3*lines['C1[uF/km]'][i], + lines['Imax[kA]'][i], + r0_ohm_per_km=lines['R0[ohm/km]'][i], + x0_ohm_per_km=lines['X0[ohm/km]'][i], + c0_nf_per_km=0) + + for i in range(0, len(nodes['Naziv'])): + pp.create_load(net, nodes['ID'][i], nodes['P [kW]'][i]/1000, nodes['Q [kVAr]'][i]/1000) + + a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type = 'balanced_positive') + # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') + +if __name__ == "__main__": + pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_2.py b/pandapower/test/harmonics/Test_Case_2.py new file mode 100644 index 0000000000..ae56030bc4 --- /dev/null +++ b/pandapower/test/harmonics/Test_Case_2.py @@ -0,0 +1,39 @@ +import pytest +import pandapower as pp +import pandapower.networks as pn +from pandapower.harmonics.balanced import balanced_thd_voltage +from pandapower.plotting import simple_plot + +def test_harmonics_case_2(): + har_a = [8.5,7.5,7,6.5,5.5,5,4.5] + har_a_angle = [170,155,140,125,110,95,80] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + net = pn.create_kerber_landnetz_kabel_1() + + net.ext_grid["s_sc_max_mva"] = 5 + net.ext_grid["rx_max"] = 3 + net.ext_grid["r0x0_max"] = 3 + net.ext_grid["x0x_max"] = 2 + + net.trafo.i0_percent = 1.25 + net.trafo["vector_group"] = 'YNyn' + net.trafo["vk0_percent"] = net.trafo["vk_percent"] + net.trafo["vkr0_percent"] = net.trafo["vkr_percent"] + net.trafo["mag0_percent"] = 100 + net.trafo["mag0_rx"] = 0 + net.trafo["si0_hv_partial"] = 100 + + net.line['r0_ohm_per_km'] = 3*net.line['r_ohm_per_km'] + net.line['x0_ohm_per_km'] = 4*net.line['x_ohm_per_km'] + net.line['c0_nf_per_km'] = 0 + + # pp.runpp_3ph(net) + + a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type = 'balanced_positive') + + #simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') + +if __name__ == "__main__": + pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_3.py b/pandapower/test/harmonics/Test_Case_3.py new file mode 100644 index 0000000000..e82a094939 --- /dev/null +++ b/pandapower/test/harmonics/Test_Case_3.py @@ -0,0 +1,46 @@ +import pytest +import pandapower as pp +from pandapower.harmonics.balanced import balanced_thd_voltage +from pandapower.plotting import simple_plot + +def test_harmonics_case_3(): + net = pp.create_empty_network() + + pp.create_bus(net, 110, index = 0) + + for i in range(1, 5): + pp.create_bus(net, 10, index = i) + + pp.create_bus(net, 0.4, index = 5) + + pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 10, rx_max = 3, r0x0_max = 3, x0x_max= 4) + + pp.create_transformer_from_parameters(net, hv_bus = 0, lv_bus = 1, sn_mva = 6, vn_hv_kv = 110, vn_lv_kv = 10, vkr_percent = 0, \ + vk_percent = 4,pfe_kw = 0, i0_percent = 0, vector_group = 'Yyn', shift_degree = 0,\ + vkr0_percent = 0.8, vk0_percent = 4, mag0_percent = 100, mag0_rx = 100, si0_hv_partial = 100) + + pp.create_line_from_parameters(net, from_bus = 1, to_bus = 2, length_km = 1, r_ohm_per_km = 0.61, x_ohm_per_km = 0.355, \ + c_nf_per_km = 0, max_i_ka = 170, r0_ohm_per_km = 0.76, x0_ohm_per_km = 1.7, c0_nf_per_km = 0) + pp.create_line_from_parameters(net, from_bus = 2, to_bus = 3, length_km = 1.5, r_ohm_per_km = 0.61, x_ohm_per_km = 0.355, \ + c_nf_per_km = 0, max_i_ka = 170, r0_ohm_per_km = 0.76, x0_ohm_per_km = 1.7, c0_nf_per_km = 0) + pp.create_line_from_parameters(net, from_bus = 3, to_bus = 4, length_km = 1.5, r_ohm_per_km = 0.61, x_ohm_per_km = 0.355, \ + c_nf_per_km = 0, max_i_ka = 170, r0_ohm_per_km = 0.76, x0_ohm_per_km = 1.7, c0_nf_per_km = 0) + + pp.create_transformer_from_parameters(net, hv_bus = 4, lv_bus = 5, sn_mva = 0.63, vn_hv_kv = 10, vn_lv_kv = 0.4, vkr_percent = 0, \ + vk_percent = 4,pfe_kw = 0, i0_percent = 0, vector_group = 'Yyn',shift_degree = 0,\ + vkr0_percent = 0, vk0_percent = 4, mag0_percent = 100, mag0_rx = 100, si0_hv_partial = 100) + + pp.create_load(net, 2, p_mw = 0.35, q_mvar = 0.115039) + pp.create_load(net, 3, p_mw = 0.375, q_mvar = 0.123257) + pp.create_load(net, 5, p_mw = 0.02, q_mvar = 0.006574) + + har_a = [11,10,9,8,7,6,5] + har_a_angle = [250,240,230,220,210,200,190] + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type = 'balanced_positive') + + # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') + +if __name__ == "__main__": + pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_4.py b/pandapower/test/harmonics/Test_Case_4.py new file mode 100644 index 0000000000..79332fa70c --- /dev/null +++ b/pandapower/test/harmonics/Test_Case_4.py @@ -0,0 +1,87 @@ +import os +import pytest +import pandas as pd +import pandapower as pp +import numpy as np +from pandapower.harmonics.unbalanced import unbalanced_thd_voltage +from pandapower.plotting import simple_plot + + +def test_harmonics_case_4(): + path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") + + # Defining the path of files containing network parameters - Modified CIGRE LV network + bus = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Bus') + line = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Line') + load = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Load') + + # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) + har_a = [8,7,6,5,4,3,2] + har_a_angle = [320,310,300,290,280,270,260] + har_b = [13,11.5,10,8.5,7,5.5,4] + har_b_angle = [300,290,280,270,280,290,300] + har_c = [10,9.25,8.5,7.75,7,6.25,5.5] + har_c_angle = [250,255,260,265,255,245,235] + + #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used + har_a_lc = [10,8,6,4,2,1.5,1] + har_a_lc_angle = [300,320,280,260,250,270,220] + har_b_lc = [15,14,13,12,11,10,9] + har_b_lc_angle = [310,297,222,200,180,144,127] + har_c_lc = [10,9,8.5,7,6.5,5,4.5] + har_c_lc_angle = [230,245,160,165,155,145,135] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + net = pp.create_empty_network() + + for i in range(0, len(bus['name'])): + coord = (bus['x'][i], bus['y'][i]) + pp.create_bus(net, bus['vn_kv'][i], name = bus['name'][i], index = bus['id'][i]-1, geodata = coord) + + pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 10, s_sc_min_mva = 10, rx_max=1, rx_min=1, r0x0_max=1, x0x_max=1) + + for i in range(0, len(load['name'])): + #We create both residential load and the LC technology, e.g., PV. For the purpose of creating and veryfing + #the tool, we determine the PV with the power of 3 kW, at the phase A at every node + #If no LC technology is used, the second asymmetric load at the same node is not needed + #Also, more than two loads can be created + pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw = load['p_a_mw'][i], q_a_mvar = load['q_a_mvar'][i], p_b_mw = load['p_b_mw'][i],\ + q_b_mvar = load['q_b_mvar'][i], p_c_mw = load['p_c_mw'][i], q_c_mvar = load['q_c_mvar'][i]) + + pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw = 0/1000, q_a_mvar = 0, p_b_mw = 0,\ + q_b_mvar = 0, p_c_mw = 0, q_c_mvar = 0) + + for i in range(0, len(line['name'])): + pp.create_line_from_parameters(net, line['from_bus'][i]-1, line['to_bus'][i]-1, line['length_km'][i],\ + line['r_ohm_per_km'][i], line['x_ohm_per_km'][i], 0, line['max_i_ka'][i],\ + r0_ohm_per_km = line['r0_ohm_per_km'][i], x0_ohm_per_km = line['x0_ohm_per_km'][i], + c0_nf_per_km = 0, name = line['name'][i]) + + pp.runpp_3ph(net) + + thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = \ + unbalanced_thd_voltage(net, harmonics, har_a, \ + har_a_angle, har_b, har_b_angle, har_c, har_c_angle,\ + har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle,\ + analysis_type = "unbalanced") + + # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') + + voltages = [] + + for i in range(0, np.shape(har_vol_a)[0]): + for j in range(0, np.shape(har_vol_a)[1]): + voltages.append(har_vol_a[i,j]) + voltages.append(har_vol_b[i,j]) + voltages.append(har_vol_c[i,j]) + + pp_volt = [] + + for i in net.res_bus_3ph.index: + pp_volt.append(net.res_bus_3ph.vm_a_pu[i]) + pp_volt.append(net.res_bus_3ph.vm_b_pu[i]) + pp_volt.append(net.res_bus_3ph.vm_c_pu[i]) + +if __name__ == "__main__": + pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_5.py b/pandapower/test/harmonics/Test_Case_5.py new file mode 100644 index 0000000000..6f2b3115f0 --- /dev/null +++ b/pandapower/test/harmonics/Test_Case_5.py @@ -0,0 +1,78 @@ +import os +import pytest +import pandas as pd +import pandapower as pp +import numpy as np +import cmath +from pandapower.harmonics.unbalanced import unbalanced_thd_voltage +from pandapower.plotting import simple_plot + + +def test_harmonics_case_5(): + path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") + + # Defining the path of files containing network parameters - Modified CIGRE LV network + bus = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Bus') + line = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Line') + load = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Load') + + # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) + har_a = [8,7,6,5,4,3,2] + har_a_angle = [320,310,300,290,280,270,260] + har_b = [13,11.5,10,8.5,7,5.5,4] + har_b_angle = [300,290,280,270,280,290,300] + har_c = [10,9.25,8.5,7.75,7,6.25,5.5] + har_c_angle = [250,255,260,265,255,245,235] + + #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used + har_a_lc = [10,8,6,4,2,1.5,1] + har_a_lc_angle = [300,320,280,260,250,270,220] + har_b_lc = [15,14,13,12,11,10,9] + har_b_lc_angle = [310,297,222,200,180,144,127] + har_c_lc = [10,9,8.5,7,6.5,5,4.5] + har_c_lc_angle = [230,245,160,165,155,145,135] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + a1 = cmath.rect(1, 2/3*cmath.pi) + a2 = cmath.rect(1, 4/3*cmath.pi) + s_base = 1000000 + matrix_A = np.matrix([[1, 1, 1], [1, a2, a1], [1, a1, a2]]) + + net = pp.create_empty_network() + + + for i in range(0, len(bus['name'])): + coord = (bus['x'][i], bus['y'][i]) + pp.create_bus(net, bus['vn_kv'][i], name = bus['name'][i], index = bus['id'][i]-1, geodata = coord) + + pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 10, s_sc_min_mva = 10, rx_max=1, rx_min=1, r0x0_max=1, x0x_max=1) + + + for i in range(0, len(load['name'])): + #We create both residential load and the LC technology, e.g., PV. For the purpose of creating and veryfing + #the tool, we determine the PV with the power of 3 kW, at the phase A at every node + #If no LC technology is used, the second asymmetric load at the same node is not needed + #Also, more than two loads can be created + pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw = load['p_a_mw'][i], q_a_mvar = load['q_a_mvar'][i], p_b_mw = load['p_b_mw'][i],\ + q_b_mvar = load['q_b_mvar'][i], p_c_mw = load['p_c_mw'][i], q_c_mvar = load['q_c_mvar'][i]) + + pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw = -10/1000, q_a_mvar = 0, p_b_mw = 0,\ + q_b_mvar = 0, p_c_mw = 0, q_c_mvar = 0) + + for i in range(0, len(line['name'])): + pp.create_line_from_parameters(net, line['from_bus'][i]-1, line['to_bus'][i]-1, line['length_km'][i],\ + line['r_ohm_per_km'][i], line['x_ohm_per_km'][i], 0, line['max_i_ka'][i],\ + r0_ohm_per_km = line['r0_ohm_per_km'][i], x0_ohm_per_km = line['x0_ohm_per_km'][i], + c0_nf_per_km = 0, name = line['name'][i]) + + pp.runpp_3ph(net) + + thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = unbalanced_thd_voltage(net, harmonics, har_a, har_a_angle, har_b, har_b_angle, har_c, har_c_angle,\ + har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle,\ + analysis_type = "unbalanced") + + # # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') + +if __name__ == "__main__": + pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_6.py b/pandapower/test/harmonics/Test_Case_6.py new file mode 100644 index 0000000000..09d114036a --- /dev/null +++ b/pandapower/test/harmonics/Test_Case_6.py @@ -0,0 +1,82 @@ +import pytest +import pandapower as pp +from pandapower.plotting import simple_plot +from pandapower.harmonics.unbalanced import unbalanced_thd_voltage, unbalanced_harmonic_current_voltage + +def test_harmonics_case_6(): + net = pp.create_empty_network() + + pp.create_bus(net, 10, index = 0) + + # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) + har_a = [8,7,6,5,4,3,2] + har_a_angle = [320,310,300,290,280,270,260] + har_b = [13,11.5,10,8.5,7,5.5,4] + har_b_angle = [300,290,280,270,280,290,300] + har_c = [10,9.25,8.5,7.75,7,6.25,5.5] + har_c_angle = [250,255,260,265,255,245,235] + + #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used + har_a_lc = [10,8,6,4,2,1.5,1] + har_a_lc_angle = [300,320,280,260,250,270,220] + har_b_lc = [15,14,13,12,11,10,9] + har_b_lc_angle = [310,297,222,200,180,144,127] + har_c_lc = [10,9,8.5,7,6.5,5,4.5] + har_c_lc_angle = [230,245,160,165,155,145,135] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + for i in range(1, 7): + pp.create_bus(net, 0.4, index = i) + + pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 5, rx_max = 3, r0x0_max = 3, x0x_max= 3) + + pp.create_transformer_from_parameters(net, hv_bus = 0, lv_bus = 1, sn_mva = 0.4, vn_hv_kv = 10, vn_lv_kv = 0.4, vkr_percent = 0, \ + vk_percent = 4,pfe_kw = 0, i0_percent = 0, vector_group = 'Yyn',shift_degree = 0,\ + vkr0_percent = 0, vk0_percent = 4, mag0_percent = 100, mag0_rx = 100, si0_hv_partial = 0.9) + + pp.create_line_from_parameters(net, from_bus = 1, to_bus = 2, length_km = 0.3, r_ohm_per_km = 0.203, x_ohm_per_km = 0.08, \ + c_nf_per_km = 0, max_i_ka = 275, r0_ohm_per_km = 0.812, x0_ohm_per_km = 0.24, c0_nf_per_km = 0) + pp.create_line_from_parameters(net, from_bus = 2, to_bus = 3, length_km = 0.1, r_ohm_per_km = 0.606, x_ohm_per_km = 0.083, \ + c_nf_per_km = 0, max_i_ka = 144, r0_ohm_per_km = 2.424, x0_ohm_per_km = 0.249, c0_nf_per_km = 0) + pp.create_line_from_parameters(net, from_bus = 2, to_bus = 4, length_km = 0.3, r_ohm_per_km = 0.203, x_ohm_per_km = 0.08, \ + c_nf_per_km = 0, max_i_ka = 275, r0_ohm_per_km = 0.812, x0_ohm_per_km = 0.24, c0_nf_per_km = 0) + pp.create_line_from_parameters(net, from_bus = 4, to_bus = 5, length_km = 0.1, r_ohm_per_km = 0.606, x_ohm_per_km = 0.083, \ + c_nf_per_km = 0, max_i_ka = 144, r0_ohm_per_km = 2.424, x0_ohm_per_km = 0.249, c0_nf_per_km = 0) + pp.create_line_from_parameters(net, from_bus = 4, to_bus = 6, length_km = 0.1, r_ohm_per_km = 0.606, x_ohm_per_km = 0.083, \ + c_nf_per_km = 0, max_i_ka = 144, r0_ohm_per_km = 2.424, x0_ohm_per_km = 0.249, c0_nf_per_km = 0) + + pp.create_asymmetric_load(net, 3, p_a_mw = 5/1000, p_b_mw = 6/1000, p_c_mw = 4/1000) + pp.create_asymmetric_load(net, 3, p_a_mw = 0, p_b_mw = 0, p_c_mw = 0) + pp.create_asymmetric_load(net, 5, p_a_mw = 6/1000, p_b_mw = 7/1000, p_c_mw = 8/1000) + pp.create_asymmetric_load(net, 5, p_a_mw = 0, p_b_mw = 0, p_c_mw = 0) + pp.create_asymmetric_load(net, 6, p_a_mw = 6/1000, p_b_mw = 8/1000, p_c_mw = 5/1000) + pp.create_asymmetric_load(net, 6, p_a_mw = 0, p_b_mw = 0, p_c_mw = 0) + + # pp.runpp_3ph(net) + + # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') + + pp_volt = [] + + for i in net.res_bus_3ph.index: + pp_volt.append(net.res_bus_3ph.vm_a_pu[i]*100) + pp_volt.append(net.res_bus_3ph.vm_b_pu[i]*100) + pp_volt.append(net.res_bus_3ph.vm_c_pu[i]*100) + + thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = \ + unbalanced_thd_voltage(net, harmonics, har_a, + har_a_angle, har_b, har_b_angle, har_c, har_c_angle, + har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle, + analysis_type = "unbalanced") + + a, b, c, d = unbalanced_harmonic_current_voltage(net, harmonics, har_a, + har_a_angle, har_b, har_b_angle, har_c, har_c_angle, + har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle, + analysis_type = "unbalanced") + + # y0 = net._ppc0['internal']['Ybus'].todense() + # y1 = net._ppc1['internal']['Ybus'].todense() + +if __name__ == "__main__": + pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/__init__.py b/pandapower/test/harmonics/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pandapower/test/test_files/harmonics/CIGRE_LV.xlsx b/pandapower/test/test_files/harmonics/CIGRE_LV.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..6ced212d2c92d7aeddd2d9a7c0d8854e9acc8e72 GIT binary patch literal 12363 zcmeHN1y>x~(nW$h1PK}(g1fuBySwY)F2Q|pf(Hxk?hxE1xQE~jgy8m(y!Y0AN2PigBp+e!TZfncmWl z0FnU=AsH)?He2G1(D}aVj)EPsK}$#0LHaDL4+>pdqGw7-CxH#R8)3h(b3K)U6`ni8 z#MMg$jV8?Oe{B@`;0~`E)+Bk$efJ6s?D-iIOz|%s+Mvcr`s2myWM3Q*;l)D@oXu(0|L2+ikKOT4Q!fX~DfBTShJqv?!iKKr*WyuyW!!}%+lW>D0;E<^8e%?@;jML1 z;-jkK2182tclbSwtgP|I91am*ZL^j~p}*!KZFH{)O}%q+g?UHeoFeX2zSWQBHh(#P zl`bLeN#)iN_pYqDBu{2|gH(L>T(kytl3oKB4xF%4QOWv?SpD5e#0%@{ zb0NXNwqH~Y<>kqE*fO{~INKOIIN1D@uyQp!hmV}7Pku9xFj1h_vK|_%Pwk# zEmtxHBj74NP?;8;^->iFMf&pFTje(!l zQED2o$u{{FQxK15#^T9jw);j8wsDWC8F85Xl^#oSA0*9|A28Xw+d0l~@xhE0amKm< zUHFVY;v0Um*7CUP0Q_BbB0WjjJcn|rEGF(dNxF|oQW;llC@DT78+PI5Y_^)<#p;s{ zKI6??L$~s|N#+IAU-93s_K27@+Z%F?pucCLBE!T?aOtOk<6gbLEX%N|k>G*}`S|9s z*eFQYD6q3z14jxSXyz+n8jGv+qk&U*Ombww_B>hR4W?q#u|sQf-09vBGhI85+*6vE zgx$1x}9RA2zQB&u3vFY7hzDdi5e>3()Ou|aO6 zi@^kjB<0oQUV&JcR=)4PYwR$jzMJWqD4C0d{nLvPd>|f9(pAQ-f`nUr zc=*<;CFypeIQ!v7o*whX>K~Fjm>P^TnFeg7IxH9bQc0-U`KPYDET#GWc>W2RXmk;iqw%4t&XSS_Sa(5{k+Dkd>`Y@-IR8e)P>iUlI&c z9M97WvyQ%Q3mdm>s!@1q`;Kg3f_eQoE5OnY148i1b9!*h+uC0_dAJf@v10Z2iq;{1 zoLbRipIfvZ+IRN2V0dxWzdQ3f+#?mli<{5AI5aL8EclBv|4a1#?9%@cK=7BP>C3(U zzmN9BN$VE@7yt!#1W)@q0qxC_xI3bSs&uZwD3rFR`q@*iC>;b6MOp=i>GQCI$sT+K zrivf8ib4I@(AHUM@#NI7^@U~HD^xd}9sIn(Z@)CoWYv`z3I>Td`8WCZ$paoli!RCX zBgvscmp}ytGlo^vPK0)ubM!wlOy1xgxO#p;8sU{|Vt$wpFfO8wO%8G_XMZan95?Z7 z)7EM)IplynVT#HXc{uf!CQX}caF`;Z|7}&@SB}u(RKr@lK2sUS6*i-x z{eXF%ur%S;a_cYpL7v)}p<(o)e(;Dfm%zuJj4ONvIEGUb>OuGAu*@~%bA*4gAZK$e zQ}dVA6d5`g7|x4n{96tHmS$!y0ES;Drk}cynWp!m3$H^#8X3<4%_xS#{Ncc8X`uC< zapAPJJBZo|UA6m7u2tFm(>31}F*^=UZ?3xCF+e{%SU^AXav_eD$)1(LIK@uZn5#&r zE}UTJxSi?w>BzH|H6|ai&Lpim+I~zyJI$cUz)z>E?xCSsQ5!cm9ux0EIb|W)4tqMJ z$jGRMW;}PV|2r$~mWf|ul#P7fM@ean71zv`tWYQO{N8u|qQP16PNLVeUNTnd`I7R9 z)vunm(G8;7v~l@JQbTW#R$S3*e!!zMIH9NJ2+y(3SYDoF*6{jB<({H!ha{XQ^u-}` z&EDIIdLPMK;qR6k7SBMK;ia$a-ZwggC8%02JjWxRA7`!`^9M+eJhG87Cu&!se|XAY z)QvnJOs-nr-WhmsuR**D+qB~MyvlBmX!_AE1-4*%ee4OyIrd?2d}3)jGC?FAD`R%{ zT250ib$#1u<*;B3@rW^doqyHUv$um_Qq{hLw>k(O5k$ed&C&G_$|2PgUHNYWA1tbgu7N zM}&zNkJJb1r*spo^^rG2=$pU_3S)h1>Z@s4kn=3 z+BsI~sXT{zOT1K!%VS=@RO6dQ@bl6dUf(I4TMPpj+;An|2Em`M5ER|*O!`tRvUP%y z%^6Q+|3VQgKITlmiT)rp4}I^0c#n@=)~GHdTsz6;HfC98a^K&8CC|1|jNDb<;L^=m zxNK+;ZDauZFZ#}@MW&DAu_t}UIK`@X#}5@dTfXMhzv#*Cl}!et*?Afjic;<`_F|o+ z8+uOTjIN>D2CxI(n=V{E{dr(n-S-=?aoZH$D~+pgqGbG5w-^)+2m#)JhAsHg?hNO; zP$zzwsY2GV09GGW9S@%U-EsmkLzvf0)h;wpGbZ(>NwdPNgP$TIQo{qvN)o2l8b=PP z^1tgyoMbgWb7p+Pe;1LxmKa^D8;&KmUo@(wiq?9fu0nTN|*ZNae%`b4VE#!WIUmsVDSI~O1IsWEsx8^q$fgR-cVq0*<*;$eJ@RpDk(wg$bY@DQM)>;R zUD-$2gEHmpM z%fDh11YAT$EUslyho)7-a-5BHCN$f7+9`&LZ(Sg!t3bz@{s~nUaxTAnTtfU09G$mU zO2_3K`tzPG28(0owo#xu=PXbT_LYHurFeqfyj-h6ahvlb_7!E-i9=J@lkTkk z#Q;Q@Z%1{L93ZXDKFiKC;fuWIrnlu%l3D`$S)-TRk;vy+*L#W26PJ7r+gMheO@)2! z%;!`bA9#19(mSwxHY5`xK{L=3UF+Y%aMvQo`YR>NQWyHDXrdrBd3W2|Z*Ar5%7(th zP2Sf^=eRkp5z$Jr9DNM(=YJqF;jnuwC?=UKt)=IPft;x1-J2cFX<^|Q=PsT-;;1zW z)H!d`UF}{Ud2UfLHrm_0Hs8r%{orB_p~2=qzJPGwWyyY}9?^t5eroRivAgR+d~}!f z&A28G1Bet0X10y!>SUZdjEA;zr1ZwzXrtzqG|Kk+s)#_)fA67#-^(uI{6;l}A%6?) z?z}ty8&7sbf>EKP#}3oR36@(gPV|Alj>b7I?b^h}4clGSGSh>%S;Z}d?;4wb5(tIY zf1Uj-5N>Ove=4m?07Kt`ks{A#5GBM${j?jxPe+WNUzf3HChqDcyXpg7`*#8?4-7>e zg?4-O7RUOE;|BYYZK6UY`5970LPt@QNw^0`+!`9)+^#Fx&2aia4Ec_8HVOkem z#vT1bu8{j3FRXD~9Gcc+?xA<6RP07~*f)sS2Z9IAtC`-gLiY{a3z=}*>LTG%_6>jr zW58>Caahq;28B!T@=5{+}JjKpH_|@8&MvcHHG&O4U!P`WJu2L zqo{sG@EGr>u8hrP!p`%Qp2_JaF(`$KgfBITQgRWj#d(IWzeU<{3a+0fGkzY_Jcx(Z zzJtX;5ML`Tkzb0tnJmqx|B0n1DweywWOggL^~nw8tar>h7Q3Hs{Y?*Ey$9ToruuLg zy^=O$yYi<&8dnOBX{9i1SMGLYcW8!k#BUx#t6K!KT5PuG=Kr5OAwZvhecb6W8 z6R#M!56|-D!!!;jHbXuy;xmMMx0_z+$ZNO1GNmmu<&!d{;KId6Vxz%jNB;Qq)=p zX}k61ti}+9dusIP&@F_2HR@ z%&_3oR@L2!fCo#ja`K$^uu2C3#CDls2!meFzCQ+n@b13S~U7S)L5OYT1f6>(Z(=2B1N98qpir%7;`ZEAG;JV zQ1G0q?PHE^I5i^3cQUCE@0ngHx@LfKL4Z2uM#~=Kd={`MCuNw%C>QQ(Y!-sHcXM#A z9Pc4B4u1}}MI`wRl3VGiZ_Tn%_KfpqDvNGRKtZfc>aL3}LnxPTdTU zBa(-8n)+5lmr2s9#W<8@GZI{RumQn|yAoRk|9D211r2|%Ip7=TnRT%)vrxZgv6lRN z7eo}AGy>M_2XM{kSng^)pV`1i?ylhz>J~01$(&mUtIJG~jm;p(i5mvuxb&%nXDm+9 zu8^&8LmQ(Rq>j{$J;F*Az8eiMOg0s}+0rIn4b&IA{q?7p%*EgJ!XF5p`fDeI_+Ryc z<(FQluf)FOID8GPp25!vgL_u<$v(OW)1%N-xVe8@BbE=sWkfi2X;Qxc=@?2@r)Li% zN|T*Y=ah0MMP>MUjk4^FlydKWUhN-ZM^YzyNW*MCn{oT5gFoidAi4#W(|^hj0#kDFgA9gUneX``N`9c=Wa zDYi0_4z;?3g%QY_d0*lA9atwU35+yT6@0iPr#s4~Qy|P6wYq#zmPgM{6?BrUQ*Cc5 zC`zE5eQepp{BPadJ1Ll1+hGbu0DOu+emE=LacGt);UG?%p#h z5{3HbN~;S9km^KL$omewQ~)KwKVhsZfz4TMHuVQ5$Fq*HEf4=dBBxg*B^kxJ=^Oer zsrMg;t|Bl7^v~vX$6~45w;^v`xUXh&)32pA||X``GX38DDcMUw{wqS)ff zc@lE^+aq3An?jnb;j%`vq8%hxXr_#h6I|u?Kig%oUT^lc-r}rZKe5|MJ!`#@aLZ>% zzlrKP)%5CwbwpIP6O)Yt*{D-0=5OM1(b>qf76s9pWBRjgx{74FN`2@HE}oX)!LBQ) z(FVE_XPU<_@5o3kSwm9YSz?kB{1-)PK|8} z+kiUN>OwXx_#}g_VNS6j5j;6m)sK1Pgu9&wba2&{O!nUXIAJ1K7X7tH5aRqw8%yB; z20A{`>ikCZbmhgb7WEVrU$hVSz(04$?N#3xQ+$KITj>q!O`I=zx9JGpNAPa7OsTebuTXj*L8TJBUs`)R^pc9K-Xz9&s9_ zty^82d+6GKj^SmLUhA8h-M7=g;M=0x<}lP^^o-Y<$mGs`R0rvgLS+E$Iq?DclX$c= zaq?nRyq=CtGZwE`iWM;QuO{z{5ymzVq2cj_boia<@50Q_TI4RHp8qk)*@KPRDGvh% zHt_amt?t*Z*2U7y&Wz#L_pc4@vDSDL9>?o0>^lKOfcpd6P8|9A=A>=n8kuQU5`IJD zv64CqGq4>G1A?6EQkkltFjd%rFKI>)5~<@y92{Bw0ns$&_#$b=F4+fLX>B*8D2bS34+(tuxJXED{x4G zH~K-wir&7=#H7ba$Q(ibpCtR#uVKB;nZAjy`4P2Md5uCX?D8{*0(F1zX!3Ma8w*zd}9>i|!DLPGQ)r zN{u9~c2mV4_*6t1*~rLTkgR7Ip^2w_g`obhECt_ z-wzV8H+(Pm6ygY$SMUka&WIyw&MLP|Kon6%fynpg-~w6wjKcjuL9~6hlGD*42z{Nm zhrO=VC%r59ngdYnyj?a+twiT7WU=X0;x$mKCedyAGcn!SC${C#}guBxxZW}eT-l+@^D9E^V9pLPdEE$NkxmAOQj*a$kK zOusTa`}astR5Fl>dPP!S+?mk6m-)sWc8M$01&_Fcd{3y?&kmGv-k)iS!eL3^b6bp2 zm(vVQ zs9k0nf|~(LZN5-Dr^3>hU7sh}Ra-bX1$kgmMWAgI+~SI^N$YX$W~zE+#UjgXvO(+; zeT|x>UVz$X1h#ca9mT+ES)_Vlo%-lG()c;l4H_D(^PP0QdgEjL=IvF})i|n>O%Sh9 zC@ef8_F7Ol=98{yR++Z#Y0`(q59}zAuf^3nFcQDUe4dP`rF?K7A#?`^KtsjpoZlyL zf1onk6Ai3`yC$Nw@v_Wb?REr~1HO2w0rcl@p6^S|XZ883%qDEoFSG|PuP%kG)7yOl zM<_tw_p=q+Yj#&}dny;SvJp~y8U2y>;u%$cgVm%s3cia_=Erjy!D{gjS)$e zF51`T1Fcug9r^96x@)CEQ!UNn7|ELBDY@b(86XS4*0*@Czm%>y?V>rYc$k`}E$*}n z&uE8<>9KpYDb}$YCbTRRYL1_@XM$N=pWYr_2|tl8?v=>$DIHwVn1SVOoQuj*tCOd- zt5s=-SOFbrv`bwK;|yaE_7*MAc*?!jKx=x~9J4};aSNU^%V7)(nvEChzfQ-9A^srP>q?ed+WY!IAXcb9y!C66idn- zrRbKQW}4m07Ew7itw2MBqp+0GlXpbVn^p>W)uFJ)K{P)CfVTQ(QS|xU;han2K^V_h zq_XcKbG7dZ?H}dooGWeUTaeM`J*C^X=6s_ddl>0*#+c=_+U5Hq>*|zN>RD)E5f$dx zZ|q#fIr+{jDFueI#2e9a2fD?q=|YlGqKS zRui>iU^m*iDCwU4sBK}^q!;OuaE-c9RfCoD&>G7^+;G%KAFFi_UEiT>8YU>|BF`S( zf$JJgiMd%-4KC21?(EWyhoD78$r_AdH*DQ)cj_+rIDydKX5BZ8XA*QL6+dS(>M$3GXrY z@hy#s;V!kWX(!mVbz9><%dcl%fNjA`5ndS@3=HXC@(XbBv^4|#Tzwa)+1acxqCW92 z`lH^Rj=1g6!ZnW0%qi{<Nz>sS z-ee%0Uh!dx*Xpi1B(b~D1z2OtSw^&4S{u)$n*i7EN#ed1POC4R!aGJ(+UVD&y6l9BIslcP<+KCBM^x@1 zPfEZP-b-btXj=y3OTuIw@kL6ri>@m7EQH=yXGqgu1@m0R>Ys>Yex33Lk!24W)i{PD zI2mJlLu?MI#7{DPgkQ`Jeyn$R+Yuy%XL#4_FY17}CgH@)EqM3dv-L|+&@lh?&i;r> zk2Q=a(uVEhp$lY}%?W-vy~$Pwyv)HuJpu6gW5y%<#*goD#E#6G_7 za-128Lp}Isp)rF*ACVRvLdEoaBT&Lgm$3^J4iXl4`f?s!F&VjFVC(H$9CaT984cCf zTA=bYx^K%CIQEjDxRq)$8hiP&_|0!@{f#jl%zPs6dfk^VQQz>dpy+?HJOorv@V%s+ zG?i>TLjBk21jaPn#qh^^cWKjf>&Vh5N7dqG?NiItm{4}kt`_#eJ^C&|8)dv9Q4Q4) zgm3_=E5+@drPQry^W2R|swqN1&fx5s33`Z-6M=Y-skB9ptyk77G4#)62d18>uZP3E zzmHgukSVe9a%TI~I_fg_%qx|B6gFW9f3?!}S}0aJ{vt={3Gbiwm9&TG5oIsyOraNZ ze*JHAwls1!GgWbMwz9YQWz{~(q*evsrIq2xeN{uB#Z@8^qHCj=@;rwls;W zeJ$tJn!k20ye5vC#5G7@fB(ca`PizaeFw8bgOMELsvYX#va)(DKRnPR4b?G8DONMo zaEv(o#f>qKB#L?%dX@nC(hD5wlg!4vSxL-X>6nP#x30bZSn=anlN;unXqF^-MQ*VJ z#6mz4z2aT?47qfUBQ{j5lRIM5prbz`$B<(fV$+ag9AX-8OibpW;Z8Z0K`_yUAycb{ zIi7aS?2Z$qZJAgVetwp{B`8a&hCL}5HIX?F|7%(j_3qn|f-x!>D+)7tE?8D-gb|@4 zDbE7@{){B@-CO}F&-~0lWW5FRT#1Daqf~nE&-jLU^xCozjfeT4X$>1HQ;a&EuAV88 z)%!Vu1a+N`o*W{yvqdT0oOP)RfFO1D8XCrq{;e!ZWCva+OGYQdN^IqQ(0Ndpfmf3K znD=>zyhNmY)%1a6!a?Qjy*8mg_~p-ThrcHc?l+1p6klGh!%MtEdnuusI+!RrJ2(Ou zOdOoee#Wbp)}jAZPG25fG*Irv!e76?mRu$N&hG@AnbDSd=QoF$91jaoQ+K@Zg`|;m z<;UJD-RR&gGCSCJ2q{-%_q-QhZ97xds-a9m)VZYdaiDV3z6n=nw+9`6cTQ;)tgn?apA=X}sO6tl3gsA@2LU&tl_?7Mbk5B<^Tajv4zl?5sp7r- zV6zCjMb={5wfc{NHGAw6aXZ8ufu5?-u;NB*D%Cm~=#n8$(lIdC<=2HNz3MJ_t5lSo z{4Bgv8GPz*9tAW2Wqon?>s<&;t4s=`%+TQiS#eWd zis2s>QfD@u5i5E_Cb7?`7G}Jp8YVz}lyFMHgsnOTWuKQ?{ol>dB?MkL`Na$cU#_D5 z*$j;w9sg&AFXs2}BQsIRVTBPf1O$0Sgt&&+GOmaiSe`G#@uS)Z9O!1DUSyKjNG2Jb zaJ&%Ls5GeMz8tenO32Trba=YRt`)|l1%h|#7!nJexS477CC+x+sN|2LMQcXK1y!R` z%5Ka1KjgAc!$g}hnb~eB6Ih~qrzus2# zljhOXw8JOn;SgvLJXTp9$W5kZyBBwNYp;b(X+ril%?Xr8@ZUz(R!IcM00Bx&YgpsN z+S~`a!?#LyNUY7~CVE1Nio*aAK9V1KcHIL7Y7%A0NVjR^_a98!(3TL$j-BAAh3YrM zYwl{!8Rz6AUtcm18K}faN=Xsod=TNv$y&Hd)JKNY)kVI7CWS>v4d;ms~zO$&f`T|ApCgsUY z#m**N^Lwyd6${%ODGbqvbyL0I6_Z=`3;lifUC#&9yF0f827rFBmcsDG;n5#uJ8%fP z7is+GR_Xt|y8oR2p>bMK=C1^Q?d1Fq@aH-I#XtVk*ZDi}_wKeopo1^s{adfw@8G|7 z)BFJi151Mc1^zF+HNW%x-c#`h)0dab^ndHH_?_kV>gOLU&S?M0@~am5JIn7m@jqDd zvHy|fSBCs|mf!Pff3TFlte9SA^m}&gcY@z{1b-0Fkp3e0V}I~F^!L@*AJ9AUe?fm= zi~Y{CWHr#^=4> z`~3y)U3-kZ#vWrobL_p=^UOKdoFD}_cwE3E01^NIpahuiXIUA-006P@000gE307ay z(ZSWi!PQ8^%gF-#oXykTj^Z6WEJF?e7W)2w<9~PsN)iW^JAs(bE@W>d*IDJ~tHh9a zz6bVVGpUNSbtUu^ndxU+TR&vKzQL4ui)X`Ii9Is!!+A7hRpnq`7Z}#lq=ps|(ATDE zNW#nB-MdHINi4r$)RblvPkzP-%Epw}{fvF!VT>T4r$jCiTTIum{@aVLMnTf z<@deWQY!qsOo4v+efUe$f6=^0q(c{UZS2;Q7TZK^z1pk zD0r*TzeFwlk89s7Xd0SBXhg96RC_=oI77u<;*@c zDsJpc&Jmej+MF|gy&CSL!7Q&XUgXZ=zQg{&2=anOXoUKwNX*@ZdTXHqVS!5MF?3}- z?btmWU2IJq9c_Qatvrndhgl%*r#Ceo6}~Dsg0YZ-x3ltNqkbBAQ;QUo352UWm0Z0r zMMW_c^|L|W`6z1=Gs>dcMEU98w1S5{eBD0{an*eLh$}Chn%O|;OjIn@ry1<8ERzMu z9S5dWD9Sxq@iAOIY43gfg5>=HzWf3uAl4#S-7I|6d4O+&pIJqtomMAVe)DsiR(WRC zXeO=2Kx02w=*O~ybw$~Eir{M=LYqch4)UUgfktlhux~_MVy%kLRO_s~ZiyQ*x0{N# zx|rAjrj(o}`?jFSV?NUIym7lY@oD{WNquxt^82okfI9ks?X3cXF1;ihYh}Q!wJdqJ zia}^>5UyPe_)s9PEcZ%pns{H=`XQ?1Ryg?qUV7dm7b{k)Gzl~QLb_%#o1{Ca`hDM} z67}=0@KwY?`}drYoLcFFQiR!3m$-E{osSK-iY^Uaqm@jYSij0@HyqX`ruFSF`m8A4 z$13CsIk>PM@eTxvSU_z1(ON;{uP{8ETeQJE=g*47R<_a~%7VX*>g2G8;Ys6%+Nm0H zB)pS1Fqegu#AX)s%Rmb!oAti6IUjI&JrPk_Sp1{5clB zSO~rl~%>D>vIWYhU2Fmc?Bf(#J z{%2HxfhK8Co&0Z)5;Y}-E+A$L>U}8CBf}jJYsQtGYEN_bF;agG;~X6&m*43U3Dd{t zI#UYlaE=hKgP~sc3m%MBIGnR~=7Ml+cyAoL13?7Gp_3svB-8JuPXc8SaIp7wc8kBG zVx)s{8b#v!Sh8{XzFjiV5{2e|pprUm9EWh*3I#_+LG_a<= zo+6#ma^s;IU#f}z}VIVCG%zbvvZ{5=37f8uVLHv+-tNA(qe@t zMLfd35Lb%OJ{fX<#(40?yR;_FgB4ean*qkwxnq{W(*w)({@%$Y=uZEWGOFlSo*_`W zgh8DmA@q$uoFdrD!on5I{@cXy!#ieb4n{-0;F!JqtPXzPKd)#cD_+gj$S_B?EFe=_D7_7oQ>UxAq$C#=moiqh!K5qNT?jd^~ zXFxh{@ZK>H`J2>9mJ_09Aj%RwW0*|IAv}Hy&YlB2cF&j#?Kxbf-U)-;0gJC%96aDLq;%@viC=9lABqag$cRZpScHcV>+*+3m$;cc~>yWh2UA$ z)deY2Po>Kp8}2vKOGvN>3oJ-BGU6*zpb6ykeJ%2+)qaf4-W`hb+pYU+W&WiF?^KW= zN(5cUx2%fS^P1#Q`1xi#b}$}&L{v!B?2!bWFM~4p#->SHd5si!Z!J)UD#vs$wY?F`C4^2~-ZqQ`d$?LhJN66f_x z`+1+Z9*wu8ddsAucEl$eeF7Z|XGTt^(tJtwx8-Bt_EnKhjKs1)cDO=ixmB2Z*U)Ck7qCPZ*(=*Y{RY36lq_g56`%RNS}7rQ{H0-_ zDahv~CnfbBdXl`qPhqpl3P+p-dL+Y=c2Onm+sd`hDv7X=Ewz0z%GqO+Y!1OEt3Klw z4A|IyR#&uln*{aj_+wghQX7wSB^`1eI2#p>eRamb$J^~T;}l3idOq)s;0g?S&vx(h zTIqN8>D>dB;{X&LUcncuPFv zh-xf)@t&Hxq7D`BhF;`8rwHH`tQw8nYJWGH$|zD*Lj(i=a@QQ~7bI*)2sWK@g9anI zyyj)ONWlFE>hL-sXI|2Cid&1Mh|y)1dB$k_lG-2yxSA>M;%?Qw^xHHtq(^6XBDgSc}yVK0{sr4Us&Nr&!82X>Yf_E6bBF4rwH?P zP8a`n8M-od8Y%B^#*=F^mSBZ%0-@yx2}_cKsyOEBk?T6{EU6k0)@{-g>BFuqhRG_MBhy;Q z30wm!A}nTw>ZWUlUULZwoD?R_HW8}3M)nZ;)Qf<1MoqW|Hyg*asf?4z_T)>L`-U~o zGL2y@;)`sj7+q43m@O`j674HqI64P<;xoL)TFSd3mRgX1&bN`N;hhzYp|s2cyFtDl z@rCg#!WJ;Tsa0+_3Q6}xiZ}kg?{g<+Kga4`W1~gk@1Nj_d~Jv*&g;*IFjZ{>g-jqw zC};!=fnj*}22xqYdUpGX8Pge0AH(5FYqj7ctVeZ^yso0Z^%!{K0Rw)75Tk!|lPHkE zV6iC`P>p;+&SL9rl|A3)6kh_K^3ntwPF_CT6j^>T6e+hDwoNp43=6Sqlg_6-=I zJ2>9TR%x!>;G#_2A4KSGlJIy5-#iwH2rQUW{Yc+v12ar@na-2SGM`Azy5>wLO3yhP zGm54sX|d3j8LoJy;mJjpW0n`{#*Fu7wDd;##4{%9i}H)cl8@QJ}VBfBu(5j!^^KpYmO{K8BUY-PT(p`g8`Ycb5X>&s?U6M_Ofxa zE^2v}$4VXPS{(%(2L=+Xum}z%-{|&yObeahRB1GAByr)|iF%A>`;3adA%KofeY8u^ z;i}dKsh)|N+x7})i?H*n=c7jI34MZh>J9Eha@h`-9wx~|1XZC@i2?=*$^(x$;!=GQ zyI74_{gcqS$4KAVGi$%X6wzBt7h@*RaK>~M?DOPor?!SH5UOS?lLxG*mpR=RbuW>X z(9Do_w2()|w7>~$@?z!b=&&;*M!$X)5#+Kn#1W8Dffu)N<*{l!VD+^*e4IZqnpz+N zsA78WaHq`bQfAB6h=D!nCEvU<@iGFg9mx9TGpCYnvvOy6b+zhT4HpX%n##n} zOM5qIKH;M>deOeD&w6fREW7eqm)8UH&6Adh9sxq$PX2uPYFi))#}sOY>v594#HYp) z^9j09NK4K9^sEnljCP0%^3$w|#3rqn8bNArvFeIL+zGV??L-I4U*`iXimaR3fgD|XcbgE)T44tG z7UkU0y(xhZ!&3bNbnk_VF;bz!7HV{i*+~b`E-y;0OT9jl4XuS=u5JNhlDuF+J;`0p z9R;(Q@do2d^G|>ao2I&le?3mSVV9IdLo-K4Xy%Cid-w&rdf8cke?LqwIc5Mc{od3( zv>yxdjK;vHeMXBtJanfTicI zK{Z8JvOZm5%Q(7vg}N3EHFnhGzq%%=D@l{|Zwp?R}aECWx7dxN3S$f4yc*%YA4*IA9x3dyv5P&SH|Acf7P<*s%<$z@qBn8GGoE7PIo zT}R4zP$oJcBjb7=H9=k33y=k2Pp45wf3QTfWz1K>C;pLBX01f{wP&;< z6wy{ju_F00hSO;N=jseOl+l=0Eo&XlXNVoIbZa~DTo+burQA}$!?3)nmFfGzJ_3%G zo)w~%;RoivB4;57wYwCYLeu0eYGV_{pY<-rBz)S#-61b7H}8%+ecg<@+F1;S1N(Ut z0==pX5c=008ujo_fEWjUQQxl7)fa>%y--f7Dqt;_NF#)kFZsGSgv>DFZOrJpx7Rmb_IpX0Z~&&OQnJ7H;( zRH!Gq2q=Ye_BkyfDVw3i`XH8Uy(M}HRC;rrkAtY4EJe52`dd;E!P3TUpm+9rm_qKR z#M~6@6T_3(vHE#R`MtAuMP+WM0b)GXG?h+pHhG5}tJmy;o$W-x9Nr)@XobdJ7Qu8ojrYwGiQw}+O`^Rv9o=oTZY*&aD?{W$ULtR{r zz7Uutt427;2=aE^ps-7)v!TmEKM?IHFAmLTvte*+V#Z)lLCBv6fRp+o=;lU?qY#V!%ym6%KreI0r*KzHUtc z4D&XpabeuEU@C-RsS2Oa#UjLASqZzw?V_eDFVEXcS+{D(ms64w1d^1=5eDX%b|nVW zc&l%3eUHS=W$83}dI2GPfN+w#`>ff#X>lsV)XHaV65LV?dh#{sHIE3IF`5g_Y|x!E z+&QOetDOkGE1(9fGmVtf?jcJcOg{7A#!XoJnE$|2`i|`Kuu{kstwfWk|8l%G2TZFH zbv^C2s4MB~b2mAI%Q2Psa#sIPemEt(u-@1xKDTq3!*Ozo_;I~&PFa);=$IRHVGIHC z_Rvlg<-0vjHXRtfTj$fJx-r3g^jYGRjvpt?8F6MJvi2Z0{{G1quBU)t`6^@L=nVtG z@`WVBXjr!4kf3(>jtFqqcpT!bmZRc3e@ynxF6z!M!Ii|h!CRZN8t?17Vt+|~!_1MJ z-C%Y$p*2bE9zn^&^yMp%JxM~&>SzoyIpR)sr-C1_(Oz*{{H(~UiE#18`FhXWtEG<} z>^Cr@*#D(|${)Hdzx!AiSXO8l`S%9-pPTpR{4ectkiy>${JoX;C-BFd0#(Rg8hgJ2 ze{FC4gf>Gnv|pMXzk>f>IsOR+0J4#P1OI>O$iLe8wVe3V(irOh`-p!O7k{<#Yh~xB zl}L<#t@He9;Md&rrvX&_-wgcBRDXs3njQRvc0uc-&}9Vwnk)Qj;qQU=CmsN3A_4&Z n5q5ut|J@<}6|P457x*9U5u|_s-CY0x1$w=LMs{-YA5Z@WR~O_j literal 0 HcmV?d00001 diff --git a/pandapower/test/test_files/harmonics/Test_1f_lines.xlsx b/pandapower/test/test_files/harmonics/Test_1f_lines.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..2cedd8b27369d6d046004738ec153711090c8cfa GIT binary patch literal 9952 zcmeHt1y>yF()Hl(9v~1Lf)j$f``{8Bh7cSE8G;43!QCyvA-Dwz5;Qmjw*(0i+zIf` zk^AKw&iQ`9y*+DnuUXx-S9kTUr)pO<2iQA^TBPQLh6P_1Z zT8c{1Lh*v`J4z}fUXQL}gN*5?52Ka3b~S1zy}I4#3x*S6(u#QoeM#|ecI?;sYILR; zMNt(lM(F_&I*(PiCCkFZ zr0CmNHcR2JW_Lvfy7b|dnB7Rt`5}QCY>yO|bz582xp{Qx-e-U*wb>IdqL>M*fz5T- zR0=zT6NI@&)z)N!QdHuv%HfLbxv4ODLeOxlxd@M`U%o-+pj~PiY~n`mo`-mw?#hw2 zWPPh_a5MQS6fhnO)y`F$uJLHLS&zXuNduqkJ?np05zaPf0NPW%K)ZvGg1pZhS>K+J z_>X?4&gP0e@@xx21&83iS^_M${J=Y<09U|&M>{&z;UhUDxk<#1p5o_^$(#qIc$ksI zRSKSABKl+QNx4|ho$N#)Phs3shMmQL zg{g@P$PUcG{@eb4;rU-2lYe>i@^}@sZcfaQ1Non!{TI_qvAB{-9uo2&o@n{LQCP&O zk1C)eU+Q3_z}0#bh#=?J=KFJSVM#b@r~k?M*JtIC_=G~V4IULCDc8y95Ozv$lEM-l_c}fE-w6c>&($%=5?7EN82;NbI5R0Vw8}zE^FPmRh z!c9wR{-_A4YUIn^Nf^!Wn@KL*#u5t?P~MwJBkBj6znrV|?zek#cJ)YG(^}wFm3g*{ zAdRPqrQOF9>5LA7TW>bi%zh1;XZV-gqjG)pxo3U`b*jp_l$9YadpPagT(<_l-6uRKz zJOeXp`Gp{4dy`m}1XA<0*|Z-3?dkjW-JmfIu;v@+f$k|P2n8Tt7Np+>qz0&HRJ zvGhcy`U={(2%7coZY4|uuBs)>P2>>cb>bkjMe3d%5mJ|Wb7sncI+kKitg!q{(%ko8 zogLx|`4t3R-1-tl#v~auT&3PL&vFSamI+7=u|^y{O5=rZ*g5(#s_04}gS8s`%|(vo!t=Uvd-I2hQ!L?8Y5`&GO0d$p$0*7s!z_a`YO^Ld_YgG z!w!pn>jawsNBljIUBRB;FiftynXDrU3#jW{Hd3OF+M6)0);|Wddt#Csu?HL>kPln& z#4U-q4J|Sw6?X{9=+18dx{p|@2pN%HtvYY}reAMrwa5BVDe%+}lUwM^cNSM0|XE18S# zLY7wH9LitUWAn(eG!zvp=W>}l`fpbpeA zFEBje@jYE3XK69g`=Z2w;27e0Fx=;U!HfMBf#9t3X>kNT;%frC10f{G;gewmRP!Gd zj|1e82=MoIcFTTXV5fr#nnmOLS+fcGH!qp!$wCVXX{1k^CqkatjBvVf+knR?2Vv{s zP~G~Q>mbd&ih#PtM(T8CtBw0+sfBaooB&U%@e?vM1z{{t!1(v6Cr>XsWzJELTYp-q zc#hb17F=U)P?RY>F69;R4soRp?U$qSV~z)J=ycuNfSU6@KIz6yjALyH2 zfxYQ}Qbs)Y>;s*9wiKeMMnMV{cC=89Eti0}8^1Uah z0)(5!y(#mkR<;7kE#-qxNQ;Dv@B}r&9g1?qRL%`s2eIA_0zwRBtqZj943g={ z)6HHT2W-BfGqlv}tk8Q=gG<=AA*FBA+#piJ{n})#d7#ZUNbhDHO@7X8x zO%A^DSgaSam3c!|kR@V-u3C5|9ahkeUOWvg?V3v#*k6}<-Nz{;l@tX|jO5Okps>zF zdXN;)@K`$Sqn$RmU0JLO9wUKZ5_G3K&M&1jMmdQADRdEgPTwcy-29UJAkeaaAT3!b z=moLd_ZTB81X4<$n$aTUo5;e1)K2u>a!iQ?)Fpwn+R*atx8*rG(0$PEJP>jwsi{a9 zOp49z>RHL}Hyy3SIrouC31v7e++65deDawjM!(h_yM=J`T0Yb)D?piXYnNfQwS4=1 zSJ@lVWACm%6S{YOL^@s?)}33p83wM z*he26EmTRotS#J0GU}_%NWAPUj7c&Yt25Zor+R;~H&|N~us7IQNRW6rSbQU}@`p0A zx#5{qyn9hPtGUx&bN=i#WYR*ng zY03)MC{5z0>hisIvL@q6$ZA@6^ z0QZ_T$6j&+A?uO@Ry#iQPj7e?pHGn}>R{E*ZzCB`a_E)ovqvQ|fEa^9=t>=4WjIwp z@@_2V#+w4zzbx9XTKIQ;X;jfm3<6ItiqlTFa%kCO1Y#glHhAC56|YL7eI}GfJsXl` zwUFbJ+C8o3zD4=B`06CVZC8bz>81WnI{zMDuGS!X5XbNP-$QFpe>jqykFbN}S_~8H z@ssyk4E^%zs9nMmon=-cMSa7bx(*L_d@DHtJiWk~2GiSjDUyyNi4)=oSZ${Ap|h!UlgWOz z{7lch)#MYQcpx-~D-713N((u@E_n6>sC}~tSJHNQ$y5pg3oOK+P($LHK3ehnuP{LS zOhunMaW@S@v zy>amIcBjn*HEN)N*EhKV9L>t_K;g$=WUyGahQW>~=M>AVlR2iil# z%VQ=+wat{JF~KvrT7Fd_RG2NRG8peAm_-e{L$j=kUhdB15}Lf|o>BT?3~W5L1xhQp zG=;DS*V}kB&l6d8eUBzs#tf8bi)S8HkHsq6)6i7bqv6{yfbMgP0^Y!D&?$T!_=Kr8 zhFa4U$t1Aavp_$IitiYL&tTf5#f+t>efeJ0|E`cWqJfk9ZIY2`xE{I2ImEnYK)K=P z)7)u3YL09Xue-zTlh6(kztitm-6l27^(;ybDD;||`_nNd_qVqs%O)LnhdW1n+f$53 zT|FH?uZB`O?vDh|s=h7~akTp!A8#j+toWR5sl`yvFHlgX9zF@LKCE1`IADl0^T)n9 zf)mT?;gsx&7svbVUc5il4{!W}dZ){+>POcCg1d_=oVZINN=8#d5&K7EiA3SB<%-D$piY1b%oSf9Hf|N3|4` zvh0IY*v5Q zhOS#HnkJ-moAgA+uxro=Wpyhwt&NJrHJ~aIG%wsRTR$AmBP4iIk~IH;RMQpOOX}4q z3f3Dl;~Cm&9@nKYO`_aWEPvWRqJ5TW3I~cWwf(~El7jjaba9kuU-i<_IRHY;6ppv@ z_Ue)JlA;m+MyA&5tY~bNRbKcFs`bc8<{6~#obewz6n3Lf^x<8$rT) z)Thz#8pfXi2{RbHij=yKCwN`KNJW1RiLi+SU(!&R_bsUPqS-ha$ z!q{vBH$ron&YQ}*m`Ka^&6z=rk$XO74AVdgwDd7ELitP!!o!ej@ix@$DN)W?#f|C- zBqnN7^<{H;OL_?gQv4N!*)U2fMX~U^fylvf#ni6vcP$V76lA8kNm4UqD{^%=Qgst- zNkZe7LB?}wdgHEoA|lAf`KnCYMySNRM2Hnd;+6$bQC(Uqv=V(J zP4;yHPe~da(42#ZI>uFVuEp8Y#_d&Uo6K7_TBvJn6z2qI0LdDw&~WmNes4=!=oGhF zvw1VQ3(rmz4xX(H4P%o(1B2#RkC4Mvy$xz33+*%8D}wK&-7^s6X49?6?qmwt z4woKg$z&uoq0)(hh6$>J$Xs!$UWq+yCTxC5SkJ~O^6j7MzQGkW_?9mIlq$m+*Hx$= z!re*xDrAXNJzJH^e@(N}>8|wC3S~Ln97R_fRa8tHg5VY(-djCAj;AQm;cp@XU3P}K z{4=VE;x?{4zM2kN&$LEN2qZ?+3Pu9m^HVLeyLiJZKaZ)Mo0Btj=A zTh3>075!G#?ugo2^@TbfR#Z&2sppsWZnFF$N0p3X{aK$4+{9UT6|*j{2NzqXtx!Dt zg1Yp#5)#Cj|b#ciDR*!eCs}8OvxFpeJ|dS5T|mp z;t@iJ>o#9_jp`wv<{lMW&5lyO@C8ahIkwG;k zRq7|cs~e9JHb>+7g5ryD)pqf}}Gs1&Obm znfMCLT_RnQ$@d5}C7#Q?5+xuhN)vGBQeJnLd~#%H)I>-n641T?jZK|gMFnaAR-BMA zS&wUdHor4(1jjz^Hmr0Wzn}RTjSYoP+?P@#F_yEvxQ+L{bUJSL`QlLMN!^!0eXX4= z7itFFGu?<$Q+H^+vzG zo5*!(?WeR`3U~ybPrW*QKgdha(F#%}RvmF*=`DH|d{Dnj%`H4j)uu5%RrcB7VqDUz zGr~RO<>l7xaksadNlz!M;Yh#$uTp?#jUm#&xuaMJ_#a#sp2Ap;2~3y~yPW6@nt~Wsle=_G@l3`lj6*kXg%7nY#Qm4`vQ!xnJ`)t z{z+VGGZ&DhrmM>2Al%!8B!WIr_Rd8$>GNHb|eE9 zE}`yze&(-A-Ip^b`K}R*kBx5r)fqze+)AA=CrHAaLHma@e(PxbDN_Gvar~01zc)C7 zVmc$bI0hFYAHY;H#VF`H=bwQ2$zss zVaKsfwZI3jTaZ{+9eT~Pbx2h6?<^9cu?feTO|rMABdDq>hWEwl5n)sZtJB?*cKi0s zRY#+h#ENdtDTI4+4qQc&kCm6&*dLpq)N!cRAF!_b=mrXkXz1VVMf*?*I}~(`$RCde zDG0=xiss7;^vnD>eq%1N>v5iK>w}G0X#aWWib@Xdg&V`H$#XZ*thPB0$0 zf?PB~uCBikg5$w|E_U-*4toZ}4wK@zcV84_wCQm1@x_XO>|yBEikHoKEEL74O!;yB z72FZQtn}_1#VTInhb)GkXvRS()?S+**}Q9y7_h6Q6cI>k2rB%(d|OI$v|`L27Yz>e zX?QNVmybvVcOY4Q^cn%@Fm+3?v>1+A_zjDsPv^{ZWX>8BQ)YS45L@>bpL;^@s5x7VhEA}Gea4&83XR>drCS8jJ&XyYuaCN9sof>)+7C+<5jdro;d9@!Oh z^O{l;`NzhEH$_6oZG~!iG8Bahe9y~h@PNS;xp>W?b4*yXy{UZM8%5guo1`bV%9~wu zVETk+{sF_X=iIX~-;u?_j=HeU7u_HE#brrtrU`#yx=cwMKK;semD}1MN5EfBN(>kI z5uV!g+p=&J>-{Nrw#+m#PUOdY+Eq^Mew2&Vrzu=(lH7d?;dsoDr+ zN7qhClbx#%m*7!-Mi!c=Rh+Kc=A_M{pN?PMof>V{>wEn({GFSk*QRsli!@evV*h5^ zs_vMP;eNIEBz?8o-EI9B`4dX-Qw06G4K1ueqpwfiSM#6oK7xg{>mn2gVo-phMDM^~pWc;@;2v>@0C z2JB72%2xK2(=<+kCSw^AqM4NS>=&zb&taEAUXaRI_bOLi&IKq0>f6SNv*RX^Ieo(xTQa2}I z$N|#5)PPS$T@sg6hW}tC1*#6=fuu_VQpPhn=cT*0SA4ra0UX^o*`bV#CB$A!3_?!E zZiTX*vOsuQSB*54sO)r*pj6&YAP+}Ea^W;uO z*PC10Y?H%?-u7#1ycKxWpGhpxM)RKNj%$=5S9IW3iFJ}2|NSfV2^V+SpNLiCmp6W< zR4@|F{Qp*`Hfxyov1sC*h|RlpAc>gQ(29L_Eu+ z3HMVo3;(7c-^SK=^|X?q`xOf3+aC6lJ^B}QjW=6Pr;~A0!|iAKIA}(g!Fud2aS~P4 zE@rx6MaR3JNaCKP0m-Nk2)_cD98R~?Y!b{-!YNF9iQzEDaoud0$DfZqGCUMKEQ~JU z?XjO9(3=_Q@~LxVLyRvI?ntg;L?zN9SsDi_NAG6e&D2G>uXb-8k8xG$faty{%>;}j z3G?nW?6%g1T%KQlsAV=UXZD(Xai}e7-E)mD&_r01n$5#TT%9qM&E4F7GqwKzK&8I}{QdiZ{}}#uEP_SGU%n=IX!vlZ z`pfh)Om{t)us$^Y`_Sr_DF9FmD~105Bdmuw4<{MFkQOlh-%b2uuJI7%VZZhliV$q% z2D_Vw-P?x%4;AMxfEME40e Date: Thu, 9 Oct 2025 11:22:36 +0200 Subject: [PATCH 2/2] refactored the code. --- pandapower/harmonics/balanced.py | 22 +- .../harmonics/harmonic_impedance_creator.py | 64 ++-- pandapower/harmonics/unbalanced.py | 206 ++++++----- pandapower/pd2ppc_zero.py | 4 +- pandapower/test/harmonics/Test_Case_1.py | 47 --- pandapower/test/harmonics/Test_Case_2.py | 39 -- pandapower/test/harmonics/Test_Case_3.py | 46 --- pandapower/test/harmonics/Test_Case_4.py | 87 ----- pandapower/test/harmonics/Test_Case_5.py | 78 ---- pandapower/test/harmonics/Test_Case_6.py | 82 ---- pandapower/test/harmonics/test_harmonics.py | 350 ++++++++++++++++++ 11 files changed, 513 insertions(+), 512 deletions(-) delete mode 100644 pandapower/test/harmonics/Test_Case_1.py delete mode 100644 pandapower/test/harmonics/Test_Case_2.py delete mode 100644 pandapower/test/harmonics/Test_Case_3.py delete mode 100644 pandapower/test/harmonics/Test_Case_4.py delete mode 100644 pandapower/test/harmonics/Test_Case_5.py delete mode 100644 pandapower/test/harmonics/Test_Case_6.py create mode 100644 pandapower/test/harmonics/test_harmonics.py diff --git a/pandapower/harmonics/balanced.py b/pandapower/harmonics/balanced.py index d4544f6ae5..d26d4f01c0 100644 --- a/pandapower/harmonics/balanced.py +++ b/pandapower/harmonics/balanced.py @@ -1,14 +1,18 @@ import math import cmath import numpy as np -import pandapower as pp +from pandapower import runpp +from pandapower.auxiliary import pandapowerNet import pandapower.harmonics.harmonic_impedance_creator as hic # formatting harmonic voltages in table and calculation of THD -def balanced_thd_voltage(net, harmonics, har, har_angle, analysis_type): - harmonics_voltage_0, harmonics_voltage = balanced_harmonic_current_voltage(net, harmonics, har, - har_angle, analysis_type) +def balanced_thd_voltage(net: pandapowerNet, + harmonics: list[int], + har: list[float], har_angle: list[float], + analysis_type: str): + + harmonics_voltage_0, harmonics_voltage = balanced_harmonic_current_voltage(net, harmonics, har, har_angle, analysis_type) # Nodes need to be sorted 0, 1, 2, 3, 4... Node with index 0 needs to be referent node (External grid is connected to this node) @@ -43,7 +47,11 @@ def balanced_thd_voltage(net, harmonics, har, har_angle, analysis_type): # calculation of harmonic voltages from harmonic currents # at the moment it is possible to define only one harmonic patter which is same for every harmonic source -def balanced_harmonic_current_voltage(net, harmonics, har, har_angle, analysis_type): +def balanced_harmonic_current_voltage(net: pandapowerNet, + harmonics: list[int], + har: list[float], har_angle: list[float], + analysis_type: str): + delta_harmonics_voltage = np.zeros([len(net.bus.name) - 1, len(har)], dtype=complex) harmonics_voltage = np.zeros([len(net.bus.name) - 1, len(har)], dtype=float) harmonic_cur_val = [] @@ -60,14 +68,14 @@ def balanced_harmonic_current_voltage(net, harmonics, har, har_angle, analysis_t z_ext = har_ext_matrix[h] if harmonics[h] == 1: - pp.runpp(net) + runpp(net) current = [] current_0 = [] s = complex(net.res_bus.p_mw[net.ext_grid.bus[0]], net.res_bus.q_mvar[net.ext_grid.bus[0]]) - current_0.append(np.conjugate(s / (math.sqrt(3) * (cmath.rect(net.res_bus.vm_pu[net.ext_grid.bus[0]], \ + current_0.append(np.conjugate(s / (math.sqrt(3) * (cmath.rect(net.res_bus.vm_pu[net.ext_grid.bus[0]], net.res_bus.va_degree[ net.ext_grid.bus[0]]))))) diff --git a/pandapower/harmonics/harmonic_impedance_creator.py b/pandapower/harmonics/harmonic_impedance_creator.py index 858c9e99ff..206f3b5d09 100644 --- a/pandapower/harmonics/harmonic_impedance_creator.py +++ b/pandapower/harmonics/harmonic_impedance_creator.py @@ -1,13 +1,15 @@ -#Script for creating impedance matrix for higher order frequencies +# Script for creating impedance matrix for higher order frequencies -import pandapower as pp import numpy as np import math import cmath +from pandapower.auxiliary import pandapowerNet from pandapower.pd2ppc import _pd2ppc from pandapower.pd2ppc_zero import _pd2ppc_zero +from pandapower.pypower.makeYbus import makeYbus -def harmonic_imp_creator(net, harmon_order, analysis_type): + +def harmonic_imp_creator(net: pandapowerNet, harmon_order: list[int], analysis_type: str): harmonics = harmon_order a1 = cmath.rect(1, 2/3*cmath.pi) @@ -18,8 +20,8 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): harmonic_matrices = [] harmonic_ext_matrices = [] - #In case of unbalanced harmonic analyses, values for zero sequence system need to be defined. - #If they are unknown, they are assumed to be same as values for positive sequence system. + # In case of unbalanced harmonic analyses, values for zero sequence system need to be defined. + # If they are unknown, they are assumed to be same as values for positive sequence system. if 'r0_ohm_per_km' not in net.line: net.line['r0_ohm_per_km'] = net.line['r_ohm_per_km'] @@ -30,7 +32,7 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): if 'c0_nf_per_km' not in net.line: net.line['c0_nf_per_km'] = net.line['c_nf_per_km'] - #Defining parameters needed for LF calculations + # Defining parameters needed for LF calculations net["_options"] = {"mode":"pf", "check_connectivity":True, "calculate_voltage_angles":True, "init_vm_pu":True, "init_va_degree":True, "consider_line_temperature":False, "voltage_depend_loads":False, @@ -41,34 +43,34 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): u_base = net.bus.vn_kv[net.ext_grid.bus[0]] - #Calculation of external grid's impedance and the network's impedance matrix for every harmonic order - #Already developed pandapower functionalities are used for matrices creation + # Calculation of external grid's impedance and the network's impedance matrix for every harmonic order + # Already developed pandapower functionalities are used for matrices creation for h in range(0, len(harmonics)): ppc_0, ppci_0 = _pd2ppc_zero(net, None) - ybus_0, yf_0, yt_0 = pp.pypower.makeYbus.makeYbus(ppci_0["baseMVA"], ppci_0["bus"], ppci_0["branch"]) + ybus_0, yf_0, yt_0 = makeYbus(ppci_0["baseMVA"], ppci_0["bus"], ppci_0["branch"]) ppc_1, ppci_1 = _pd2ppc(net, 1) - ybus_1, yf_1, yt_1 = pp.pypower.makeYbus.makeYbus(ppci_1["baseMVA"], ppci_1["bus"], ppci_1["branch"]) + ybus_1, yf_1, yt_1 = makeYbus(ppci_1["baseMVA"], ppci_1["bus"], ppci_1["branch"]) ppc_2, ppci_2 = _pd2ppc(net, 2) - ybus_2, yf_2, yt_2 = pp.pypower.makeYbus.makeYbus(ppci_2["baseMVA"], ppci_2["bus"], ppci_2["branch"]) + ybus_2, yf_2, yt_2 = makeYbus(ppci_2["baseMVA"], ppci_2["bus"], ppci_2["branch"]) - #Full sequence admittance matrices + # Full sequence admittance matrices zero_full = ybus_0.todense() pos_full = ybus_1.todense() neg_full = ybus_2.todense() - #Removing referent node related row and column of each matrix - #It is assumed that first row and column are referent node related + # Removing referent node related row and column of each matrix + # It is assumed that first row and column are referent node related zero = np.delete(np.delete(zero_full, 0, 0), 0, 1) pos = np.delete(np.delete(pos_full, 0, 0), 0, 1) neg = np.delete(np.delete(neg_full, 0, 0), 0, 1) - #External grid impedance calculation + # External grid impedance calculation z_pos_ohm = u_base**2/net.ext_grid.s_sc_max_mva delta_pos = math.atan(1/net.ext_grid.rx_max) @@ -82,24 +84,24 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): x_zer = z_zer_ohm*math.sin(delta_zer)*harmonics[h] z_zer = complex(r_zer, x_zer) - #Depenent on the analysis type, impedances are created + # Depenent on the analysis type, impedances are created if analysis_type == 'unbalanced': phase_mat_y = np.zeros([3*(np.shape(zero)[0]), 3*(np.shape(zero)[1])], dtype = complex) z_ext_0 = np.zeros([3,3], dtype = complex) - #Z012 matrix is created for the referent node + # Z012 matrix is created for the referent node z_ext_0[0, 0] = z_zer/(u_base**2/3) z_ext_0[1, 1] = z_pos/(u_base**2/3) z_ext_0[2, 2] = z_pos/(u_base**2/3) - #Transformation from the sequence to the phase system + # Transformation from the sequence to the phase system z_ext_abc = np.matmul(np.matmul(np.linalg.inv(matrix_A), z_ext_0), matrix_A) for i in range(0, np.shape(zero)[0]): for j in range(0, np.shape(zero)[1]): - #Y012 is created for every other element of Y0,1,2 matrices + # Y012 is created for every other element of Y0,1,2 matrices y_012 = np.zeros([3,3], dtype = complex) z = zero[i, j] @@ -110,16 +112,16 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): y_012[1,1] = p y_012[2,2] = n - #Transformation from the sequence to the phase system + # Transformation from the sequence to the phase system y_abc = np.matmul(np.matmul(np.linalg.inv(matrix_A), y_012), matrix_A) for row in range (0, 3): for col in range(0, 3): phase_mat_y[3*i+row, 3*j+col] = y_abc[row, col] - #Z matrix as an inverse of the Y matrix + # Z matrix as an inverse of the Y matrix phase_mat_z = np.linalg.inv(phase_mat_y)*3 - #Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X + # Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X for r in range(0, np.shape(phase_mat_z)[0]): for c in range(0, np.shape(phase_mat_z)[1]): aux = phase_mat_z[r,c] @@ -128,16 +130,16 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): imp = complex(res, reac) phase_mat_z[r,c] = imp - #Appending the referent node and network impedances + # Appending the referent node and network impedances harmonic_matrices.append(phase_mat_z) harmonic_ext_matrices.append(z_ext_abc) elif analysis_type == 'balanced_positive': - #Only positive sequence system is observed + # Only positive sequence system is observed z_pos_net = np.linalg.inv(pos) - #Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X + # Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X for r in range(0, np.shape(z_pos_net)[0]): for c in range(0, np.shape(z_pos_net)[1]): aux = z_pos_net[r,c] @@ -146,19 +148,19 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): imp = complex(res, reac) z_pos_net[r,c] = imp - #Appending the referent node and network impedances + # Appending the referent node and network impedances harmonic_matrices.append(z_pos_net) harmonic_ext_matrices.append(z_pos/(u_base**2)) elif analysis_type == 'balanced_all': - #Balanced network but all sequence systems are used - #Dependent on the harmonic order, zero, positive or negative sequence impedance are considered + # Balanced network but all sequence systems are used + # Dependent on the harmonic order, zero, positive or negative sequence impedance are considered z_pos_net = np.linalg.inv(pos) z_neg_net = np.linalg.inv(neg) z_zero_net = np.linalg.inv(zero) - #Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X + # Dependent on the harmonic order, reactive part of the impedance (reactance) is calculated as h*X for r in range(0, np.shape(z_pos_net)[0]): for c in range(0, np.shape(z_pos_net)[1]): aux_p = z_pos_net[r,c] @@ -179,7 +181,7 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): imp_z = complex(res_z, reac_z) z_zero_net[r,c] = imp_z - #Appending the referent node and network impedances + # Appending the referent node and network impedances if harmonics[h] % 3 == 0: harmonic_matrices.append(z_zero_net) harmonic_ext_matrices.append(z_zer/(u_base**2)) @@ -190,5 +192,5 @@ def harmonic_imp_creator(net, harmon_order, analysis_type): harmonic_matrices.append(z_neg_net) harmonic_ext_matrices.append(z_pos/(u_base**2)) - return(harmonic_matrices, harmonic_ext_matrices) + return harmonic_matrices, harmonic_ext_matrices \ No newline at end of file diff --git a/pandapower/harmonics/unbalanced.py b/pandapower/harmonics/unbalanced.py index c2a99ca423..9ddef24249 100644 --- a/pandapower/harmonics/unbalanced.py +++ b/pandapower/harmonics/unbalanced.py @@ -1,7 +1,8 @@ import math import cmath import numpy as np -import pandapower as pp +from pandapower import runpp_3ph +from pandapower.auxiliary import pandapowerNet import pandapower.harmonics.harmonic_impedance_creator as hic #import pandapower.harmonics.unbalanced.voltages_currents_calculation as har_vol_calc @@ -9,15 +10,27 @@ # formatting harmonic voltages in table and calculation of THD # Nodes need to be sorted 0, 1, 2, 3, 4... Node with index 0 needs to be referent node (External grid is connected to this node) -def unbalanced_thd_voltage(network, harmonics, har_a, har_a_angle, har_b, har_b_angle, har_c, har_c_angle, - har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle, - analysis_type): - harmonics_voltage_0, harmonics_voltage, harmonics_current_0, harmonics_current = \ - unbalanced_harmonic_current_voltage(network, harmonics, - har_a, har_a_angle, har_b, har_b_angle, har_c, har_c_angle, - har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, - har_c_lc_angle, - analysis_type) +def unbalanced_thd_voltage(net: pandapowerNet, + harmonics: list[int], + har_a: list[float], har_a_angle: list[float], + har_b: list[float], har_b_angle: list[float], + har_c: list[float], har_c_angle: list[float], + har_a_lc: list[float], har_a_lc_angle: list[float], + har_b_lc: list[float], har_b_lc_angle: list[float], + har_c_lc: list[float], har_c_lc_angle: list[float], + analysis_type: str): + + harmonics_voltage_0, harmonics_voltage, harmonics_current_0, harmonics_current = unbalanced_harmonic_current_voltage( + net, + harmonics, + har_a, har_a_angle, + har_b, har_b_angle, + har_c, har_c_angle, + har_a_lc, har_a_lc_angle, + har_b_lc, har_b_lc_angle, + har_c_lc, har_c_lc_angle, + analysis_type + ) thd_a = [] thd_b = [] @@ -42,9 +55,9 @@ def unbalanced_thd_voltage(network, harmonics, har_a, har_a_angle, har_b, har_b_ thd_c.append(sum_thd) for i in range(0, len(thd_a)): - thd_a[i] = math.sqrt(thd_a[i]) / (network.res_bus_3ph.vm_a_pu[i + 1]) * 100 - thd_b[i] = math.sqrt(thd_b[i]) / (network.res_bus_3ph.vm_b_pu[i + 1]) * 100 - thd_c[i] = math.sqrt(thd_c[i]) / (network.res_bus_3ph.vm_c_pu[i + 1]) * 100 + thd_a[i] = math.sqrt(thd_a[i]) / (net.res_bus_3ph.vm_a_pu[i + 1]) * 100 + thd_b[i] = math.sqrt(thd_b[i]) / (net.res_bus_3ph.vm_b_pu[i + 1]) * 100 + thd_c[i] = math.sqrt(thd_c[i]) / (net.res_bus_3ph.vm_c_pu[i + 1]) * 100 sum_thd_a_0 = 0 sum_thd_b_0 = 0 @@ -55,42 +68,49 @@ def unbalanced_thd_voltage(network, harmonics, har_a, har_a_angle, har_b, har_b_ sum_thd_b_0 += harmonics_voltage_0[1, i] ** 2 sum_thd_c_0 += harmonics_voltage_0[2, i] ** 2 - thd_a.insert(0, math.sqrt(sum_thd_a_0) / network.res_bus_3ph.vm_a_pu[0] * 100) - thd_b.insert(0, math.sqrt(sum_thd_b_0) / network.res_bus_3ph.vm_b_pu[0] * 100) - thd_c.insert(0, math.sqrt(sum_thd_c_0) / network.res_bus_3ph.vm_c_pu[0] * 100) + thd_a.insert(0, math.sqrt(sum_thd_a_0) / net.res_bus_3ph.vm_a_pu[0] * 100) + thd_b.insert(0, math.sqrt(sum_thd_b_0) / net.res_bus_3ph.vm_b_pu[0] * 100) + thd_c.insert(0, math.sqrt(sum_thd_c_0) / net.res_bus_3ph.vm_c_pu[0] * 100) - harmonics_voltage_a = np.zeros([int(len(network.bus.name) * 3 / 3), len(har_a)], dtype=float) - harmonics_voltage_b = np.zeros([int(len(network.bus.name) * 3 / 3), len(har_a)], dtype=float) - harmonics_voltage_c = np.zeros([int(len(network.bus.name) * 3 / 3), len(har_a)], dtype=float) + harmonics_voltage_a = np.zeros([int(len(net.bus.name) * 3 / 3), len(har_a)], dtype=float) + harmonics_voltage_b = np.zeros([int(len(net.bus.name) * 3 / 3), len(har_a)], dtype=float) + harmonics_voltage_c = np.zeros([int(len(net.bus.name) * 3 / 3), len(har_a)], dtype=float) # Harmonic voltage (percentage) of each phase and each harmoni for i in range(0, np.shape(harmonics_voltage_a)[0] - 1): for j in range(0, np.shape(harmonics_voltage)[1]): - harmonics_voltage_a[i + 1, j] = harmonics_voltage[3 * i, j] * 100 * (1 / network.res_bus_3ph.vm_a_pu[i + 1]) + harmonics_voltage_a[i + 1, j] = harmonics_voltage[3 * i, j] * 100 * (1 / net.res_bus_3ph.vm_a_pu[i + 1]) harmonics_voltage_b[i + 1, j] = harmonics_voltage[3 * i + 1, j] * 100 * ( - 1 / network.res_bus_3ph.vm_b_pu[i + 1]) + 1 / net.res_bus_3ph.vm_b_pu[i + 1]) harmonics_voltage_c[i + 1, j] = harmonics_voltage[3 * i + 2, j] * 100 * ( - 1 / network.res_bus_3ph.vm_c_pu[i + 1]) + 1 / net.res_bus_3ph.vm_c_pu[i + 1]) for i in range(0, np.shape(harmonics_voltage_0)[1]): - harmonics_voltage_a[0, i] = harmonics_voltage_0[0, i] * 100 * (1 / network.res_bus_3ph.vm_a_pu[0]) - harmonics_voltage_b[0, i] = harmonics_voltage_0[1, i] * 100 * (1 / network.res_bus_3ph.vm_b_pu[0]) - harmonics_voltage_c[0, i] = harmonics_voltage_0[2, i] * 100 * (1 / network.res_bus_3ph.vm_c_pu[0]) + harmonics_voltage_a[0, i] = harmonics_voltage_0[0, i] * 100 * (1 / net.res_bus_3ph.vm_a_pu[0]) + harmonics_voltage_b[0, i] = harmonics_voltage_0[1, i] * 100 * (1 / net.res_bus_3ph.vm_b_pu[0]) + harmonics_voltage_c[0, i] = harmonics_voltage_0[2, i] * 100 * (1 / net.res_bus_3ph.vm_c_pu[0]) return thd_a, thd_b, thd_c, harmonics_voltage_a, harmonics_voltage_b, harmonics_voltage_c -def unbalanced_harmonic_current_voltage(network, harmonics, har_a, har_a_angle, har_b, har_b_angle, har_c, har_c_angle, \ - har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle, \ - analysis_type): - delta_harmonics_voltage = np.zeros([len(network.bus.name) * 3 - 3, len(har_a)], dtype=complex) - harmonics_voltage = np.zeros([len(network.bus.name) * 3 - 3, len(har_a)], dtype=float) - harmonics_current = np.zeros([len(network.bus.name) * 3 - 3, len(har_a)], dtype=float) +def unbalanced_harmonic_current_voltage(net: pandapowerNet, + harmonics: list[int], + har_a: list[float], har_a_angle: list[float], + har_b: list[float], har_b_angle: list[float], + har_c: list[float], har_c_angle: list[float], + har_a_lc: list[float], har_a_lc_angle: list[float], + har_b_lc: list[float], har_b_lc_angle: list[float], + har_c_lc: list[float], har_c_lc_angle: list[float], + analysis_type: str): + + delta_harmonics_voltage = np.zeros([len(net.bus.name) * 3 - 3, len(har_a)], dtype=complex) + harmonics_voltage = np.zeros([len(net.bus.name) * 3 - 3, len(har_a)], dtype=float) + harmonics_current = np.zeros([len(net.bus.name) * 3 - 3, len(har_a)], dtype=float) u_harmonics_0 = np.zeros([3, len(har_a)], dtype=complex) i_harmonics_0 = np.zeros([3, len(har_a)], dtype=float) - har_matrices, har_ext_matrix = hic.harmonic_imp_creator(network, harmonics, analysis_type) + har_matrices, har_ext_matrix = hic.harmonic_imp_creator(net, harmonics, analysis_type) - s_base = network.sn_mva * 1e6 + s_base = net.sn_mva * 1e6 for h in range(0, len(harmonics)): phase_mat_z = har_matrices[h] @@ -98,7 +118,7 @@ def unbalanced_harmonic_current_voltage(network, harmonics, har_a, har_a_angle, if harmonics[h] == 1: - pp.runpp_3ph(network) + runpp_3ph(net) current = [] current_res = [] @@ -106,72 +126,72 @@ def unbalanced_harmonic_current_voltage(network, harmonics, har_a, har_a_angle, current_0 = [] # p.u. values - s_a = complex(network.res_bus_3ph.p_a_mw[network.ext_grid.bus[0]] * 1000000 / s_base, - network.res_bus_3ph.q_a_mvar[network.ext_grid.bus[0]] * 1000000 / s_base) - s_b = complex(network.res_bus_3ph.p_b_mw[network.ext_grid.bus[0]] * 1000000 / s_base, - network.res_bus_3ph.q_b_mvar[network.ext_grid.bus[0]] * 1000000 / s_base) - s_c = complex(network.res_bus_3ph.p_c_mw[network.ext_grid.bus[0]] * 1000000 / s_base, - network.res_bus_3ph.q_c_mvar[network.ext_grid.bus[0]] * 1000000 / s_base) - - current_0.append(np.conjugate(s_a / (cmath.rect(network.res_bus_3ph.vm_a_pu[network.ext_grid.bus[0]], - network.res_bus_3ph.va_a_degree[ - network.ext_grid.bus[0]] * cmath.pi / 180)))) - current_0.append(np.conjugate(s_b / (cmath.rect(network.res_bus_3ph.vm_b_pu[network.ext_grid.bus[0]], - network.res_bus_3ph.va_b_degree[ - network.ext_grid.bus[0]] * cmath.pi / 180)))) - current_0.append(np.conjugate(s_c / (cmath.rect(network.res_bus_3ph.vm_c_pu[network.ext_grid.bus[0]], - network.res_bus_3ph.va_c_degree[ - network.ext_grid.bus[0]] * cmath.pi / 180)))) - - for i in network.res_bus_3ph.index: + s_a = complex(net.res_bus_3ph.p_a_mw[net.ext_grid.bus[0]] * 1000000 / s_base, + net.res_bus_3ph.q_a_mvar[net.ext_grid.bus[0]] * 1000000 / s_base) + s_b = complex(net.res_bus_3ph.p_b_mw[net.ext_grid.bus[0]] * 1000000 / s_base, + net.res_bus_3ph.q_b_mvar[net.ext_grid.bus[0]] * 1000000 / s_base) + s_c = complex(net.res_bus_3ph.p_c_mw[net.ext_grid.bus[0]] * 1000000 / s_base, + net.res_bus_3ph.q_c_mvar[net.ext_grid.bus[0]] * 1000000 / s_base) + + current_0.append(np.conjugate(s_a / (cmath.rect(net.res_bus_3ph.vm_a_pu[net.ext_grid.bus[0]], + net.res_bus_3ph.va_a_degree[ + net.ext_grid.bus[0]] * cmath.pi / 180)))) + current_0.append(np.conjugate(s_b / (cmath.rect(net.res_bus_3ph.vm_b_pu[net.ext_grid.bus[0]], + net.res_bus_3ph.va_b_degree[ + net.ext_grid.bus[0]] * cmath.pi / 180)))) + current_0.append(np.conjugate(s_c / (cmath.rect(net.res_bus_3ph.vm_c_pu[net.ext_grid.bus[0]], + net.res_bus_3ph.va_c_degree[ + net.ext_grid.bus[0]] * cmath.pi / 180)))) + + for i in net.res_bus_3ph.index: connected = 0 - if i != network.ext_grid.bus[0]: - for j in range(0, len(network.asymmetric_load.index), 2): + if i != net.ext_grid.bus[0]: + for j in range(0, len(net.asymmetric_load.index), 2): # The assumption is that two harmonic sources are connected to the node # It can be meodified, additional sources can be added, or number of sources at every node can be reduced - if network.asymmetric_load.bus[j] == i: + if net.asymmetric_load.bus[j] == i: connected = 1 # an asymmetric load is connected to the observed node ##Residential - s_a_res = complex(network.asymmetric_load.p_a_mw[j] * 1000000 / s_base, - network.asymmetric_load.q_a_mvar[j] * 1000000 / s_base) - s_b_res = complex(network.asymmetric_load.p_b_mw[j] * 1000000 / s_base, - network.asymmetric_load.q_b_mvar[j] * 1000000 / s_base) - s_c_res = complex(network.asymmetric_load.p_c_mw[j] * 1000000 / s_base, - network.asymmetric_load.q_c_mvar[j] * 1000000 / s_base) + s_a_res = complex(net.asymmetric_load.p_a_mw[j] * 1000000 / s_base, + net.asymmetric_load.q_a_mvar[j] * 1000000 / s_base) + s_b_res = complex(net.asymmetric_load.p_b_mw[j] * 1000000 / s_base, + net.asymmetric_load.q_b_mvar[j] * 1000000 / s_base) + s_c_res = complex(net.asymmetric_load.p_c_mw[j] * 1000000 / s_base, + net.asymmetric_load.q_c_mvar[j] * 1000000 / s_base) # PV - s_a_pv = complex(network.asymmetric_load.p_a_mw[j + 1] * 1000000 / s_base, - network.asymmetric_load.q_a_mvar[j + 1] * 1000000 / s_base) - s_b_pv = complex(network.asymmetric_load.p_b_mw[j + 1] * 1000000 / s_base, - network.asymmetric_load.q_b_mvar[j + 1] * 1000000 / s_base) - s_c_pv = complex(network.asymmetric_load.p_c_mw[j + 1] * 1000000 / s_base, - network.asymmetric_load.q_c_mvar[j + 1] * 1000000 / s_base) + s_a_pv = complex(net.asymmetric_load.p_a_mw[j + 1] * 1000000 / s_base, + net.asymmetric_load.q_a_mvar[j + 1] * 1000000 / s_base) + s_b_pv = complex(net.asymmetric_load.p_b_mw[j + 1] * 1000000 / s_base, + net.asymmetric_load.q_b_mvar[j + 1] * 1000000 / s_base) + s_c_pv = complex(net.asymmetric_load.p_c_mw[j + 1] * 1000000 / s_base, + net.asymmetric_load.q_c_mvar[j + 1] * 1000000 / s_base) break if connected == 1: # If the node is a node where a load/generator is connected we calculate the currents # Else currents at that node are 0 - current_res.append(np.conjugate(s_a_res / (cmath.rect(network.res_bus_3ph.vm_a_pu[i], - network.res_bus_3ph.va_a_degree[ + current_res.append(np.conjugate(s_a_res / (cmath.rect(net.res_bus_3ph.vm_a_pu[i], + net.res_bus_3ph.va_a_degree[ i] * cmath.pi / 180)))) - current_res.append(np.conjugate(s_b_res / (cmath.rect(network.res_bus_3ph.vm_b_pu[i], - network.res_bus_3ph.va_b_degree[ + current_res.append(np.conjugate(s_b_res / (cmath.rect(net.res_bus_3ph.vm_b_pu[i], + net.res_bus_3ph.va_b_degree[ i] * cmath.pi / 180)))) - current_res.append(np.conjugate(s_c_res / (cmath.rect(network.res_bus_3ph.vm_c_pu[i], - network.res_bus_3ph.va_c_degree[ + current_res.append(np.conjugate(s_c_res / (cmath.rect(net.res_bus_3ph.vm_c_pu[i], + net.res_bus_3ph.va_c_degree[ i] * cmath.pi / 180)))) - current_pv.append(np.conjugate(s_a_pv / (cmath.rect(network.res_bus_3ph.vm_a_pu[i], - network.res_bus_3ph.va_a_degree[ + current_pv.append(np.conjugate(s_a_pv / (cmath.rect(net.res_bus_3ph.vm_a_pu[i], + net.res_bus_3ph.va_a_degree[ i] * cmath.pi / 180)))) - current_pv.append(np.conjugate(s_b_pv / (cmath.rect(network.res_bus_3ph.vm_b_pu[i], - network.res_bus_3ph.va_b_degree[ + current_pv.append(np.conjugate(s_b_pv / (cmath.rect(net.res_bus_3ph.vm_b_pu[i], + net.res_bus_3ph.va_b_degree[ i] * cmath.pi / 180)))) - current_pv.append(np.conjugate(s_c_pv / (cmath.rect(network.res_bus_3ph.vm_c_pu[i], - network.res_bus_3ph.va_c_degree[ + current_pv.append(np.conjugate(s_c_pv / (cmath.rect(net.res_bus_3ph.vm_c_pu[i], + net.res_bus_3ph.va_c_degree[ i] * cmath.pi / 180)))) else: @@ -184,21 +204,21 @@ def unbalanced_harmonic_current_voltage(network, harmonics, har_a, har_a_angle, current_pv.append(0 + 0j) # p.u. values of current in each node. Must be equal to s_res + s_pv - s_a = complex(network.res_bus_3ph.p_a_mw[i] * 1000000 / s_base, - network.res_bus_3ph.q_a_mvar[i] * 1000000 / s_base) - s_b = complex(network.res_bus_3ph.p_b_mw[i] * 1000000 / s_base, - network.res_bus_3ph.q_b_mvar[i] * 1000000 / s_base) - s_c = complex(network.res_bus_3ph.p_c_mw[i] * 1000000 / s_base, - network.res_bus_3ph.q_c_mvar[i] * 1000000 / s_base) - - current.append(np.conjugate(s_a / (cmath.rect(network.res_bus_3ph.vm_a_pu[i], - network.res_bus_3ph.va_a_degree[ + s_a = complex(net.res_bus_3ph.p_a_mw[i] * 1000000 / s_base, + net.res_bus_3ph.q_a_mvar[i] * 1000000 / s_base) + s_b = complex(net.res_bus_3ph.p_b_mw[i] * 1000000 / s_base, + net.res_bus_3ph.q_b_mvar[i] * 1000000 / s_base) + s_c = complex(net.res_bus_3ph.p_c_mw[i] * 1000000 / s_base, + net.res_bus_3ph.q_c_mvar[i] * 1000000 / s_base) + + current.append(np.conjugate(s_a / (cmath.rect(net.res_bus_3ph.vm_a_pu[i], + net.res_bus_3ph.va_a_degree[ i] * cmath.pi / 180)))) - current.append(np.conjugate(s_b / (cmath.rect(network.res_bus_3ph.vm_b_pu[i], - network.res_bus_3ph.va_b_degree[ + current.append(np.conjugate(s_b / (cmath.rect(net.res_bus_3ph.vm_b_pu[i], + net.res_bus_3ph.va_b_degree[ i] * cmath.pi / 180)))) - current.append(np.conjugate(s_c / (cmath.rect(network.res_bus_3ph.vm_c_pu[i], - network.res_bus_3ph.va_c_degree[ + current.append(np.conjugate(s_c / (cmath.rect(net.res_bus_3ph.vm_c_pu[i], + net.res_bus_3ph.va_c_degree[ i] * cmath.pi / 180)))) else: @@ -208,7 +228,7 @@ def unbalanced_harmonic_current_voltage(network, harmonics, har_a, har_a_angle, harmonics_current_pv = [] harmonics_angle_pv = [] - for i in range(0, len(network.bus.name)): + for i in range(0, len(net.bus.name)): # Residential household devices harmonics harmonics_current_res.append(har_a[h - 1]) harmonics_current_res.append(har_b[h - 1]) diff --git a/pandapower/pd2ppc_zero.py b/pandapower/pd2ppc_zero.py index 06343ad351..af4b3bbcdb 100644 --- a/pandapower/pd2ppc_zero.py +++ b/pandapower/pd2ppc_zero.py @@ -210,8 +210,8 @@ def _add_trafo_sc_impedance_zero(net, ppc, trafo_df=None, k_st=None): tap_hv = np.square(vn_trafo_hv / vn_bus_hv) * net.sn_mva if mode == 'pf_3ph': if vector_group.lower() not in ["ynyn", "dyn", "yzn"]: - raise NotImplementedError("Calculation of 3-phase power flow is only implemented for the transformer " - "vector groups 'YNyn', 'Dyn', 'Yzn'") + raise NotImplementedError(f"Calculation of 3-phase power flow is only implemented for the transformer " + f"vector groups 'YNyn', 'Dyn', 'Yzn' and not for {vector_group.lower()}") # ============================================================================= # Changing base from transformer base to Network base to get Zpu(Net) # Zbase = (kV).squared/S_mva diff --git a/pandapower/test/harmonics/Test_Case_1.py b/pandapower/test/harmonics/Test_Case_1.py deleted file mode 100644 index efd9e5951a..0000000000 --- a/pandapower/test/harmonics/Test_Case_1.py +++ /dev/null @@ -1,47 +0,0 @@ -import os -import pytest -import pandas as pd -import pandapower as pp -from pandapower.harmonics.balanced import balanced_thd_voltage -from pandapower.plotting.simple_plot import simple_plot - -def test_harmonics_case_1(): - path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") - - nodes = pd.read_excel(os.path.join(path, 'Test_1f_Loads.xlsx')) - lines = pd.read_excel(os.path.join(path,'Test_1f_lines.xlsx')) - - har_a = [11,10,9,8,7,6,5] - har_a_angle = [250,240,230,220,210,200,190] - - harmonics = [1, 3, 5, 7, 9, 11, 13, 15] - # harmonics = [1, 3] - - net = pp.create_empty_network() - - for i in range(0, len(nodes['Naziv'])): - pp.create_bus(net, 10, nodes['Naziv'][i], int(nodes['ID'][i])) - - pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 10, s_sc_min_mva = 10, rx_max=5, rx_min=0, r0x0_max=5, x0x_max=3) - - for i in range (0, len(lines['Length'])): - pp.create_line_from_parameters(net, - lines['Start Node'][i], - lines['End Node'][i], - lines['Length'][i], - lines['R1[ohm/km]'][i], - lines['X1[ohm/km]'][i], - 1e3*lines['C1[uF/km]'][i], - lines['Imax[kA]'][i], - r0_ohm_per_km=lines['R0[ohm/km]'][i], - x0_ohm_per_km=lines['X0[ohm/km]'][i], - c0_nf_per_km=0) - - for i in range(0, len(nodes['Naziv'])): - pp.create_load(net, nodes['ID'][i], nodes['P [kW]'][i]/1000, nodes['Q [kVAr]'][i]/1000) - - a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type = 'balanced_positive') - # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') - -if __name__ == "__main__": - pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_2.py b/pandapower/test/harmonics/Test_Case_2.py deleted file mode 100644 index ae56030bc4..0000000000 --- a/pandapower/test/harmonics/Test_Case_2.py +++ /dev/null @@ -1,39 +0,0 @@ -import pytest -import pandapower as pp -import pandapower.networks as pn -from pandapower.harmonics.balanced import balanced_thd_voltage -from pandapower.plotting import simple_plot - -def test_harmonics_case_2(): - har_a = [8.5,7.5,7,6.5,5.5,5,4.5] - har_a_angle = [170,155,140,125,110,95,80] - - harmonics = [1, 3, 5, 7, 9, 11, 13, 15] - - net = pn.create_kerber_landnetz_kabel_1() - - net.ext_grid["s_sc_max_mva"] = 5 - net.ext_grid["rx_max"] = 3 - net.ext_grid["r0x0_max"] = 3 - net.ext_grid["x0x_max"] = 2 - - net.trafo.i0_percent = 1.25 - net.trafo["vector_group"] = 'YNyn' - net.trafo["vk0_percent"] = net.trafo["vk_percent"] - net.trafo["vkr0_percent"] = net.trafo["vkr_percent"] - net.trafo["mag0_percent"] = 100 - net.trafo["mag0_rx"] = 0 - net.trafo["si0_hv_partial"] = 100 - - net.line['r0_ohm_per_km'] = 3*net.line['r_ohm_per_km'] - net.line['x0_ohm_per_km'] = 4*net.line['x_ohm_per_km'] - net.line['c0_nf_per_km'] = 0 - - # pp.runpp_3ph(net) - - a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type = 'balanced_positive') - - #simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') - -if __name__ == "__main__": - pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_3.py b/pandapower/test/harmonics/Test_Case_3.py deleted file mode 100644 index e82a094939..0000000000 --- a/pandapower/test/harmonics/Test_Case_3.py +++ /dev/null @@ -1,46 +0,0 @@ -import pytest -import pandapower as pp -from pandapower.harmonics.balanced import balanced_thd_voltage -from pandapower.plotting import simple_plot - -def test_harmonics_case_3(): - net = pp.create_empty_network() - - pp.create_bus(net, 110, index = 0) - - for i in range(1, 5): - pp.create_bus(net, 10, index = i) - - pp.create_bus(net, 0.4, index = 5) - - pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 10, rx_max = 3, r0x0_max = 3, x0x_max= 4) - - pp.create_transformer_from_parameters(net, hv_bus = 0, lv_bus = 1, sn_mva = 6, vn_hv_kv = 110, vn_lv_kv = 10, vkr_percent = 0, \ - vk_percent = 4,pfe_kw = 0, i0_percent = 0, vector_group = 'Yyn', shift_degree = 0,\ - vkr0_percent = 0.8, vk0_percent = 4, mag0_percent = 100, mag0_rx = 100, si0_hv_partial = 100) - - pp.create_line_from_parameters(net, from_bus = 1, to_bus = 2, length_km = 1, r_ohm_per_km = 0.61, x_ohm_per_km = 0.355, \ - c_nf_per_km = 0, max_i_ka = 170, r0_ohm_per_km = 0.76, x0_ohm_per_km = 1.7, c0_nf_per_km = 0) - pp.create_line_from_parameters(net, from_bus = 2, to_bus = 3, length_km = 1.5, r_ohm_per_km = 0.61, x_ohm_per_km = 0.355, \ - c_nf_per_km = 0, max_i_ka = 170, r0_ohm_per_km = 0.76, x0_ohm_per_km = 1.7, c0_nf_per_km = 0) - pp.create_line_from_parameters(net, from_bus = 3, to_bus = 4, length_km = 1.5, r_ohm_per_km = 0.61, x_ohm_per_km = 0.355, \ - c_nf_per_km = 0, max_i_ka = 170, r0_ohm_per_km = 0.76, x0_ohm_per_km = 1.7, c0_nf_per_km = 0) - - pp.create_transformer_from_parameters(net, hv_bus = 4, lv_bus = 5, sn_mva = 0.63, vn_hv_kv = 10, vn_lv_kv = 0.4, vkr_percent = 0, \ - vk_percent = 4,pfe_kw = 0, i0_percent = 0, vector_group = 'Yyn',shift_degree = 0,\ - vkr0_percent = 0, vk0_percent = 4, mag0_percent = 100, mag0_rx = 100, si0_hv_partial = 100) - - pp.create_load(net, 2, p_mw = 0.35, q_mvar = 0.115039) - pp.create_load(net, 3, p_mw = 0.375, q_mvar = 0.123257) - pp.create_load(net, 5, p_mw = 0.02, q_mvar = 0.006574) - - har_a = [11,10,9,8,7,6,5] - har_a_angle = [250,240,230,220,210,200,190] - harmonics = [1, 3, 5, 7, 9, 11, 13, 15] - - a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type = 'balanced_positive') - - # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') - -if __name__ == "__main__": - pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_4.py b/pandapower/test/harmonics/Test_Case_4.py deleted file mode 100644 index 79332fa70c..0000000000 --- a/pandapower/test/harmonics/Test_Case_4.py +++ /dev/null @@ -1,87 +0,0 @@ -import os -import pytest -import pandas as pd -import pandapower as pp -import numpy as np -from pandapower.harmonics.unbalanced import unbalanced_thd_voltage -from pandapower.plotting import simple_plot - - -def test_harmonics_case_4(): - path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") - - # Defining the path of files containing network parameters - Modified CIGRE LV network - bus = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Bus') - line = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Line') - load = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Load') - - # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) - har_a = [8,7,6,5,4,3,2] - har_a_angle = [320,310,300,290,280,270,260] - har_b = [13,11.5,10,8.5,7,5.5,4] - har_b_angle = [300,290,280,270,280,290,300] - har_c = [10,9.25,8.5,7.75,7,6.25,5.5] - har_c_angle = [250,255,260,265,255,245,235] - - #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used - har_a_lc = [10,8,6,4,2,1.5,1] - har_a_lc_angle = [300,320,280,260,250,270,220] - har_b_lc = [15,14,13,12,11,10,9] - har_b_lc_angle = [310,297,222,200,180,144,127] - har_c_lc = [10,9,8.5,7,6.5,5,4.5] - har_c_lc_angle = [230,245,160,165,155,145,135] - - harmonics = [1, 3, 5, 7, 9, 11, 13, 15] - - net = pp.create_empty_network() - - for i in range(0, len(bus['name'])): - coord = (bus['x'][i], bus['y'][i]) - pp.create_bus(net, bus['vn_kv'][i], name = bus['name'][i], index = bus['id'][i]-1, geodata = coord) - - pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 10, s_sc_min_mva = 10, rx_max=1, rx_min=1, r0x0_max=1, x0x_max=1) - - for i in range(0, len(load['name'])): - #We create both residential load and the LC technology, e.g., PV. For the purpose of creating and veryfing - #the tool, we determine the PV with the power of 3 kW, at the phase A at every node - #If no LC technology is used, the second asymmetric load at the same node is not needed - #Also, more than two loads can be created - pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw = load['p_a_mw'][i], q_a_mvar = load['q_a_mvar'][i], p_b_mw = load['p_b_mw'][i],\ - q_b_mvar = load['q_b_mvar'][i], p_c_mw = load['p_c_mw'][i], q_c_mvar = load['q_c_mvar'][i]) - - pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw = 0/1000, q_a_mvar = 0, p_b_mw = 0,\ - q_b_mvar = 0, p_c_mw = 0, q_c_mvar = 0) - - for i in range(0, len(line['name'])): - pp.create_line_from_parameters(net, line['from_bus'][i]-1, line['to_bus'][i]-1, line['length_km'][i],\ - line['r_ohm_per_km'][i], line['x_ohm_per_km'][i], 0, line['max_i_ka'][i],\ - r0_ohm_per_km = line['r0_ohm_per_km'][i], x0_ohm_per_km = line['x0_ohm_per_km'][i], - c0_nf_per_km = 0, name = line['name'][i]) - - pp.runpp_3ph(net) - - thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = \ - unbalanced_thd_voltage(net, harmonics, har_a, \ - har_a_angle, har_b, har_b_angle, har_c, har_c_angle,\ - har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle,\ - analysis_type = "unbalanced") - - # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') - - voltages = [] - - for i in range(0, np.shape(har_vol_a)[0]): - for j in range(0, np.shape(har_vol_a)[1]): - voltages.append(har_vol_a[i,j]) - voltages.append(har_vol_b[i,j]) - voltages.append(har_vol_c[i,j]) - - pp_volt = [] - - for i in net.res_bus_3ph.index: - pp_volt.append(net.res_bus_3ph.vm_a_pu[i]) - pp_volt.append(net.res_bus_3ph.vm_b_pu[i]) - pp_volt.append(net.res_bus_3ph.vm_c_pu[i]) - -if __name__ == "__main__": - pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_5.py b/pandapower/test/harmonics/Test_Case_5.py deleted file mode 100644 index 6f2b3115f0..0000000000 --- a/pandapower/test/harmonics/Test_Case_5.py +++ /dev/null @@ -1,78 +0,0 @@ -import os -import pytest -import pandas as pd -import pandapower as pp -import numpy as np -import cmath -from pandapower.harmonics.unbalanced import unbalanced_thd_voltage -from pandapower.plotting import simple_plot - - -def test_harmonics_case_5(): - path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") - - # Defining the path of files containing network parameters - Modified CIGRE LV network - bus = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Bus') - line = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Line') - load = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name = 'Load') - - # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) - har_a = [8,7,6,5,4,3,2] - har_a_angle = [320,310,300,290,280,270,260] - har_b = [13,11.5,10,8.5,7,5.5,4] - har_b_angle = [300,290,280,270,280,290,300] - har_c = [10,9.25,8.5,7.75,7,6.25,5.5] - har_c_angle = [250,255,260,265,255,245,235] - - #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used - har_a_lc = [10,8,6,4,2,1.5,1] - har_a_lc_angle = [300,320,280,260,250,270,220] - har_b_lc = [15,14,13,12,11,10,9] - har_b_lc_angle = [310,297,222,200,180,144,127] - har_c_lc = [10,9,8.5,7,6.5,5,4.5] - har_c_lc_angle = [230,245,160,165,155,145,135] - - harmonics = [1, 3, 5, 7, 9, 11, 13, 15] - - a1 = cmath.rect(1, 2/3*cmath.pi) - a2 = cmath.rect(1, 4/3*cmath.pi) - s_base = 1000000 - matrix_A = np.matrix([[1, 1, 1], [1, a2, a1], [1, a1, a2]]) - - net = pp.create_empty_network() - - - for i in range(0, len(bus['name'])): - coord = (bus['x'][i], bus['y'][i]) - pp.create_bus(net, bus['vn_kv'][i], name = bus['name'][i], index = bus['id'][i]-1, geodata = coord) - - pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 10, s_sc_min_mva = 10, rx_max=1, rx_min=1, r0x0_max=1, x0x_max=1) - - - for i in range(0, len(load['name'])): - #We create both residential load and the LC technology, e.g., PV. For the purpose of creating and veryfing - #the tool, we determine the PV with the power of 3 kW, at the phase A at every node - #If no LC technology is used, the second asymmetric load at the same node is not needed - #Also, more than two loads can be created - pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw = load['p_a_mw'][i], q_a_mvar = load['q_a_mvar'][i], p_b_mw = load['p_b_mw'][i],\ - q_b_mvar = load['q_b_mvar'][i], p_c_mw = load['p_c_mw'][i], q_c_mvar = load['q_c_mvar'][i]) - - pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw = -10/1000, q_a_mvar = 0, p_b_mw = 0,\ - q_b_mvar = 0, p_c_mw = 0, q_c_mvar = 0) - - for i in range(0, len(line['name'])): - pp.create_line_from_parameters(net, line['from_bus'][i]-1, line['to_bus'][i]-1, line['length_km'][i],\ - line['r_ohm_per_km'][i], line['x_ohm_per_km'][i], 0, line['max_i_ka'][i],\ - r0_ohm_per_km = line['r0_ohm_per_km'][i], x0_ohm_per_km = line['x0_ohm_per_km'][i], - c0_nf_per_km = 0, name = line['name'][i]) - - pp.runpp_3ph(net) - - thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = unbalanced_thd_voltage(net, harmonics, har_a, har_a_angle, har_b, har_b_angle, har_c, har_c_angle,\ - har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle,\ - analysis_type = "unbalanced") - - # # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') - -if __name__ == "__main__": - pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/Test_Case_6.py b/pandapower/test/harmonics/Test_Case_6.py deleted file mode 100644 index 09d114036a..0000000000 --- a/pandapower/test/harmonics/Test_Case_6.py +++ /dev/null @@ -1,82 +0,0 @@ -import pytest -import pandapower as pp -from pandapower.plotting import simple_plot -from pandapower.harmonics.unbalanced import unbalanced_thd_voltage, unbalanced_harmonic_current_voltage - -def test_harmonics_case_6(): - net = pp.create_empty_network() - - pp.create_bus(net, 10, index = 0) - - # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) - har_a = [8,7,6,5,4,3,2] - har_a_angle = [320,310,300,290,280,270,260] - har_b = [13,11.5,10,8.5,7,5.5,4] - har_b_angle = [300,290,280,270,280,290,300] - har_c = [10,9.25,8.5,7.75,7,6.25,5.5] - har_c_angle = [250,255,260,265,255,245,235] - - #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used - har_a_lc = [10,8,6,4,2,1.5,1] - har_a_lc_angle = [300,320,280,260,250,270,220] - har_b_lc = [15,14,13,12,11,10,9] - har_b_lc_angle = [310,297,222,200,180,144,127] - har_c_lc = [10,9,8.5,7,6.5,5,4.5] - har_c_lc_angle = [230,245,160,165,155,145,135] - - harmonics = [1, 3, 5, 7, 9, 11, 13, 15] - - for i in range(1, 7): - pp.create_bus(net, 0.4, index = i) - - pp.create_ext_grid(net, bus = 0, s_sc_max_mva = 5, rx_max = 3, r0x0_max = 3, x0x_max= 3) - - pp.create_transformer_from_parameters(net, hv_bus = 0, lv_bus = 1, sn_mva = 0.4, vn_hv_kv = 10, vn_lv_kv = 0.4, vkr_percent = 0, \ - vk_percent = 4,pfe_kw = 0, i0_percent = 0, vector_group = 'Yyn',shift_degree = 0,\ - vkr0_percent = 0, vk0_percent = 4, mag0_percent = 100, mag0_rx = 100, si0_hv_partial = 0.9) - - pp.create_line_from_parameters(net, from_bus = 1, to_bus = 2, length_km = 0.3, r_ohm_per_km = 0.203, x_ohm_per_km = 0.08, \ - c_nf_per_km = 0, max_i_ka = 275, r0_ohm_per_km = 0.812, x0_ohm_per_km = 0.24, c0_nf_per_km = 0) - pp.create_line_from_parameters(net, from_bus = 2, to_bus = 3, length_km = 0.1, r_ohm_per_km = 0.606, x_ohm_per_km = 0.083, \ - c_nf_per_km = 0, max_i_ka = 144, r0_ohm_per_km = 2.424, x0_ohm_per_km = 0.249, c0_nf_per_km = 0) - pp.create_line_from_parameters(net, from_bus = 2, to_bus = 4, length_km = 0.3, r_ohm_per_km = 0.203, x_ohm_per_km = 0.08, \ - c_nf_per_km = 0, max_i_ka = 275, r0_ohm_per_km = 0.812, x0_ohm_per_km = 0.24, c0_nf_per_km = 0) - pp.create_line_from_parameters(net, from_bus = 4, to_bus = 5, length_km = 0.1, r_ohm_per_km = 0.606, x_ohm_per_km = 0.083, \ - c_nf_per_km = 0, max_i_ka = 144, r0_ohm_per_km = 2.424, x0_ohm_per_km = 0.249, c0_nf_per_km = 0) - pp.create_line_from_parameters(net, from_bus = 4, to_bus = 6, length_km = 0.1, r_ohm_per_km = 0.606, x_ohm_per_km = 0.083, \ - c_nf_per_km = 0, max_i_ka = 144, r0_ohm_per_km = 2.424, x0_ohm_per_km = 0.249, c0_nf_per_km = 0) - - pp.create_asymmetric_load(net, 3, p_a_mw = 5/1000, p_b_mw = 6/1000, p_c_mw = 4/1000) - pp.create_asymmetric_load(net, 3, p_a_mw = 0, p_b_mw = 0, p_c_mw = 0) - pp.create_asymmetric_load(net, 5, p_a_mw = 6/1000, p_b_mw = 7/1000, p_c_mw = 8/1000) - pp.create_asymmetric_load(net, 5, p_a_mw = 0, p_b_mw = 0, p_c_mw = 0) - pp.create_asymmetric_load(net, 6, p_a_mw = 6/1000, p_b_mw = 8/1000, p_c_mw = 5/1000) - pp.create_asymmetric_load(net, 6, p_a_mw = 0, p_b_mw = 0, p_c_mw = 0) - - # pp.runpp_3ph(net) - - # simple_plot(net, line_width = 2.5, trafo_size=2.0, bus_size = 1.5, ext_grid_size = 1.5, line_color = 'black') - - pp_volt = [] - - for i in net.res_bus_3ph.index: - pp_volt.append(net.res_bus_3ph.vm_a_pu[i]*100) - pp_volt.append(net.res_bus_3ph.vm_b_pu[i]*100) - pp_volt.append(net.res_bus_3ph.vm_c_pu[i]*100) - - thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = \ - unbalanced_thd_voltage(net, harmonics, har_a, - har_a_angle, har_b, har_b_angle, har_c, har_c_angle, - har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle, - analysis_type = "unbalanced") - - a, b, c, d = unbalanced_harmonic_current_voltage(net, harmonics, har_a, - har_a_angle, har_b, har_b_angle, har_c, har_c_angle, - har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, har_c_lc, har_c_lc_angle, - analysis_type = "unbalanced") - - # y0 = net._ppc0['internal']['Ybus'].todense() - # y1 = net._ppc1['internal']['Ybus'].todense() - -if __name__ == "__main__": - pytest.main([__file__, "-xs"]) diff --git a/pandapower/test/harmonics/test_harmonics.py b/pandapower/test/harmonics/test_harmonics.py new file mode 100644 index 0000000000..abe3f84ee7 --- /dev/null +++ b/pandapower/test/harmonics/test_harmonics.py @@ -0,0 +1,350 @@ +import os +import cmath +import pytest +import numpy as np +import pandas as pd +import pandapower as pp +import pandapower.networks as nw +from pandapower.harmonics.balanced import balanced_thd_voltage +from pandapower.harmonics.unbalanced import unbalanced_thd_voltage, unbalanced_harmonic_current_voltage +from pandapower.plotting.simple_plot import simple_plot + +def test_harmonics_case_1(): + path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") + + nodes = pd.read_excel(os.path.join(path, 'Test_1f_Loads.xlsx')) + lines = pd.read_excel(os.path.join(path,'Test_1f_lines.xlsx')) + + har_a = [11, 10, 9, 8, 7, 6, 5] + har_a_angle = [250, 240, 230, 220, 210, 200, 190] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + # harmonics=[1, 3] + + net = pp.create_empty_network() + + for i in range(0, len(nodes['Naziv'])): + pp.create_bus(net, 10, nodes['Naziv'][i], int(nodes['ID'][i])) + + pp.create_ext_grid(net, bus=0, s_sc_max_mva=10, s_sc_min_mva=10, rx_max=5, rx_min=0, r0x0_max=5, x0x_max=3) + + for i in range (0, len(lines['Length'])): + pp.create_line_from_parameters(net, + from_bus=lines['Start Node'][i], + to_bus=lines['End Node'][i], + length_km=lines['Length'][i], + r_ohm_per_km=lines['R1[ohm/km]'][i], + x_ohm_per_km=lines['X1[ohm/km]'][i], + c_nf_per_km=1e3*lines['C1[uF/km]'][i], + max_i_ka=lines['Imax[kA]'][i], + r0_ohm_per_km=lines['R0[ohm/km]'][i], + x0_ohm_per_km=lines['X0[ohm/km]'][i], + c0_nf_per_km=0) + + for i in range(0, len(nodes['Naziv'])): + pp.create_load(net, nodes['ID'][i], nodes['P [kW]'][i]/1000, nodes['Q [kVAr]'][i]/1000) + + a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type='balanced_positive') + # simple_plot(net, line_width=2.5, trafo_size=2.0, bus_size=1.5, ext_grid_size=1.5, line_color='black') + + +def test_harmonics_case_2(): + har_a = [8.5, 7.5, 7, 6.5, 5.5, 5, 4.5] + har_a_angle = [170, 155, 140, 125, 110, 95, 80] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + net = nw.create_kerber_landnetz_kabel_1() + + net.ext_grid["s_sc_max_mva"]=5 + net.ext_grid["rx_max"]=3 + net.ext_grid["r0x0_max"]=3 + net.ext_grid["x0x_max"]=2 + + net.trafo.i0_percent=1.25 + net.trafo["vector_group"]='YNyn' + net.trafo["vk0_percent"]=net.trafo["vk_percent"] + net.trafo["vkr0_percent"]=net.trafo["vkr_percent"] + net.trafo["mag0_percent"]=100 + net.trafo["mag0_rx"]=0 + net.trafo["si0_hv_partial"]=100 + + net.line['r0_ohm_per_km']=3*net.line['r_ohm_per_km'] + net.line['x0_ohm_per_km']=4*net.line['x_ohm_per_km'] + net.line['c0_nf_per_km']=0 + + # pp.runpp_3ph(net) + + a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type='balanced_positive') + + #simple_plot(net, line_width=2.5, trafo_size=2.0, bus_size=1.5, ext_grid_size=1.5, line_color='black') + + +def test_harmonics_case_3(): + net=pp.create_empty_network() + + pp.create_bus(net, 110, index=0) + + for i in range(1, 5): + pp.create_bus(net, 10, index=i) + + pp.create_bus(net, 0.4, index=5) + + pp.create_ext_grid(net, bus=0, s_sc_max_mva=10, rx_max=3, r0x0_max=3, x0x_max= 4) + + pp.create_transformer_from_parameters(net, hv_bus=0, lv_bus=1, sn_mva=6, vn_hv_kv=110, vn_lv_kv=10, + vkr_percent=0, vk_percent=4,pfe_kw=0, i0_percent=0, vector_group='Yyn', + shift_degree=0, vkr0_percent=0.8, vk0_percent=4, mag0_percent=100, + mag0_rx=100, si0_hv_partial=100) + + pp.create_line_from_parameters(net, from_bus=1, to_bus=2, length_km=1, r_ohm_per_km=0.61, + x_ohm_per_km=0.355, c_nf_per_km=0, max_i_ka=170, r0_ohm_per_km=0.76, + x0_ohm_per_km=1.7, c0_nf_per_km=0) + pp.create_line_from_parameters(net, from_bus=2, to_bus=3, length_km=1.5, r_ohm_per_km=0.61, + x_ohm_per_km=0.355, c_nf_per_km=0, max_i_ka=170, r0_ohm_per_km=0.76, + x0_ohm_per_km=1.7, c0_nf_per_km=0) + pp.create_line_from_parameters(net, from_bus=3, to_bus=4, length_km=1.5, r_ohm_per_km=0.61, + x_ohm_per_km=0.355, c_nf_per_km=0, max_i_ka=170, r0_ohm_per_km=0.76, + x0_ohm_per_km=1.7, c0_nf_per_km=0) + + pp.create_transformer_from_parameters(net, hv_bus=4, lv_bus=5, sn_mva=0.63, vn_hv_kv=10, vn_lv_kv=0.4, + vkr_percent=0, vk_percent=4,pfe_kw=0, i0_percent=0, vector_group='Yyn', + shift_degree=0, vkr0_percent=0, vk0_percent=4, mag0_percent=100, + mag0_rx=100, si0_hv_partial=100) + + pp.create_load(net, 2, p_mw=0.35, q_mvar=0.115039) + pp.create_load(net, 3, p_mw=0.375, q_mvar=0.123257) + pp.create_load(net, 5, p_mw=0.02, q_mvar=0.006574) + + har_a = [11, 10, 9, 8, 7, 6, 5] + har_a_angle = [250, 240, 230, 220, 210, 200, 190] + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + a, b = balanced_thd_voltage(net, harmonics, har_a, har_a_angle, analysis_type='balanced_positive') + + # simple_plot(net, line_width=2.5, trafo_size=2.0, bus_size=1.5, ext_grid_size=1.5, line_color='black') + + +def test_harmonics_case_4(): + path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") + + # Defining the path of files containing network parameters - Modified CIGRE LV network + bus = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name='Bus') + line = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name='Line') + load = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name='Load') + + # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) + har_a = [8, 7, 6, 5, 4, 3, 2] + har_a_angle = [320, 310, 300, 290, 280, 270, 260] + har_b = [13, 11.5, 10, 8.5, 7, 5.5, 4] + har_b_angle = [300, 290, 280, 270, 280, 290, 300] + har_c = [10, 9.25, 8.5, 7.75, 7, 6.25, 5.5] + har_c_angle = [250, 255, 260, 265, 255, 245, 235] + + #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used + har_a_lc = [10, 8, 6, 4, 2, 1.5, 1] + har_a_lc_angle = [300, 320, 280, 260, 250, 270, 220] + har_b_lc = [15, 14, 13, 12, 11, 10, 9] + har_b_lc_angle = [310, 297, 222, 200, 180, 144, 127] + har_c_lc = [10, 9, 8.5, 7, 6.5, 5, 4.5] + har_c_lc_angle = [230, 245, 160, 165, 155, 145, 135] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + net = pp.create_empty_network() + + for i in range(0, len(bus['name'])): + coord=(bus['x'][i], bus['y'][i]) + pp.create_bus(net, vn_kv=bus['vn_kv'][i], name=bus['name'][i], index=bus['id'][i]-1, geodata=coord) + + pp.create_ext_grid(net, bus=0, s_sc_max_mva=10, s_sc_min_mva=10, rx_max=1, rx_min=1, r0x0_max=1, x0x_max=1) + + for i in range(0, len(load['name'])): + #We create both residential load and the LC technology, e.g., PV. For the purpose of creating and veryfing + #the tool, we determine the PV with the power of 3 kW, at the phase A at every node + #If no LC technology is used, the second asymmetric load at the same node is not needed + #Also, more than two loads can be created + pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw=load['p_a_mw'][i], q_a_mvar=load['q_a_mvar'][i], p_b_mw=load['p_b_mw'][i],\ + q_b_mvar=load['q_b_mvar'][i], p_c_mw=load['p_c_mw'][i], q_c_mvar=load['q_c_mvar'][i]) + + pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw=0/1000, q_a_mvar=0, p_b_mw=0,\ + q_b_mvar=0, p_c_mw=0, q_c_mvar=0) + + for i in range(0, len(line['name'])): + pp.create_line_from_parameters(net, line['from_bus'][i]-1, line['to_bus'][i]-1, line['length_km'][i],\ + line['r_ohm_per_km'][i], line['x_ohm_per_km'][i], 0, line['max_i_ka'][i],\ + r0_ohm_per_km=line['r0_ohm_per_km'][i], x0_ohm_per_km=line['x0_ohm_per_km'][i], + c0_nf_per_km=0, name=line['name'][i]) + + pp.runpp_3ph(net) + + thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = unbalanced_thd_voltage(net, harmonics, har_a, har_a_angle, + har_b, har_b_angle, har_c, har_c_angle, + har_a_lc, har_a_lc_angle, har_b_lc, + har_b_lc_angle, har_c_lc, har_c_lc_angle, + analysis_type="unbalanced") + + # simple_plot(net, line_width=2.5, trafo_size=2.0, bus_size=1.5, ext_grid_size=1.5, line_color='black') + + voltages = [] + + for i in range(0, np.shape(har_vol_a)[0]): + for j in range(0, np.shape(har_vol_a)[1]): + voltages.append(har_vol_a[i,j]) + voltages.append(har_vol_b[i,j]) + voltages.append(har_vol_c[i,j]) + + pp_volt=[] + + for i in net.res_bus_3ph.index: + pp_volt.append(net.res_bus_3ph.vm_a_pu[i]) + pp_volt.append(net.res_bus_3ph.vm_b_pu[i]) + pp_volt.append(net.res_bus_3ph.vm_c_pu[i]) + + +def test_harmonics_case_5(): + path = os.path.join(pp.pp_dir, "test", "test_files", "harmonics") + + # Defining the path of files containing network parameters - Modified CIGRE LV network + bus = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name='Bus') + line = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name='Line') + load = pd.read_excel(os.path.join(path, 'CIGRE_LV.xlsx'), sheet_name='Load') + + # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) + har_a = [8, 7, 6, 5, 4, 3, 2] + har_a_angle = [320, 310, 300, 290, 280, 270, 260] + har_b = [13, 11.5, 10, 8.5, 7, 5.5, 4] + har_b_angle = [300, 290, 280, 270, 280, 290, 300] + har_c = [10, 9.25, 8.5, 7.75, 7, 6.25, 5.5] + har_c_angle = [250, 255, 260, 265, 255, 245, 235] + + #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used + har_a_lc = [10, 8, 6, 4, 2, 1.5, 1] + har_a_lc_angle = [300, 320, 280, 260, 250, 270, 220] + har_b_lc = [15, 14, 13, 12, 11, 10, 9] + har_b_lc_angle = [310, 297, 222, 200, 180, 144, 127] + har_c_lc = [10, 9, 8.5, 7, 6.5, 5, 4.5] + har_c_lc_angle = [230, 245, 160, 165, 155, 145, 135] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + a1 = cmath.rect(1, 2/3*cmath.pi) + a2 = cmath.rect(1, 4/3*cmath.pi) + s_base = 1000000 + matrix_A = np.matrix([[1, 1, 1], [1, a2, a1], [1, a1, a2]]) + + net = pp.create_empty_network() + + + for i in range(0, len(bus['name'])): + coord = (bus['x'][i], bus['y'][i]) + pp.create_bus(net, bus['vn_kv'][i], name=bus['name'][i], index=bus['id'][i]-1, geodata=coord) + + pp.create_ext_grid(net, bus=0, s_sc_max_mva=10, s_sc_min_mva=10, rx_max=1, rx_min=1, r0x0_max=1, x0x_max=1) + + + for i in range(0, len(load['name'])): + #We create both residential load and the LC technology, e.g., PV. For the purpose of creating and veryfing + #the tool, we determine the PV with the power of 3 kW, at the phase A at every node + #If no LC technology is used, the second asymmetric load at the same node is not needed + #Also, more than two loads can be created + pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw=load['p_a_mw'][i], q_a_mvar=load['q_a_mvar'][i], p_b_mw=load['p_b_mw'][i], + q_b_mvar=load['q_b_mvar'][i], p_c_mw=load['p_c_mw'][i], q_c_mvar=load['q_c_mvar'][i]) + + pp.create_asymmetric_load(net, load['bus'][i]-1, p_a_mw=-10/1000, q_a_mvar=0, p_b_mw=0,\ + q_b_mvar=0, p_c_mw=0, q_c_mvar=0) + + for i in range(0, len(line['name'])): + pp.create_line_from_parameters(net, line['from_bus'][i]-1, line['to_bus'][i]-1, line['length_km'][i], + line['r_ohm_per_km'][i], line['x_ohm_per_km'][i], 0, line['max_i_ka'][i], + r0_ohm_per_km=line['r0_ohm_per_km'][i], x0_ohm_per_km=line['x0_ohm_per_km'][i], + c0_nf_per_km=0, name=line['name'][i]) + + pp.runpp_3ph(net) + + thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = unbalanced_thd_voltage(net, harmonics, har_a, har_a_angle, + har_b, har_b_angle, har_c, har_c_angle, + har_a_lc, har_a_lc_angle, har_b_lc, + har_b_lc_angle, har_c_lc, har_c_lc_angle, + analysis_type="unbalanced") + + # # simple_plot(net, line_width=2.5, trafo_size=2.0, bus_size=1.5, ext_grid_size=1.5, line_color='black') + + +def test_harmonics_case_6(): + net = pp.create_empty_network() + + pp.create_bus(net, 10, index=0) + + # 3rd, 5th, 7th, 9th, 11th, 13th, 15th harmonic currents (percentage of fundamental harmonic and angles) + har_a = [8, 7, 6, 5, 4, 3, 2] + har_a_angle = [320, 310, 300, 290, 280, 270, 260] + har_b = [13, 11.5, 10, 8.5, 7, 5.5, 4] + har_b_angle = [300, 290, 280, 270, 280, 290, 300] + har_c = [10, 9.25, 8.5, 7.75, 7, 6.25, 5.5] + har_c_angle = [250, 255, 260, 265, 255, 245, 235] + + #LC technology harmonic pattern, e.g., PV - does not need to be defined, or more than one LC technology can be used + har_a_lc = [10, 8, 6, 4, 2, 1.5, 1] + har_a_lc_angle = [300, 320, 280, 260, 250, 270, 220] + har_b_lc = [15, 14, 13, 12, 11, 10, 9] + har_b_lc_angle = [310, 297, 222, 200, 180, 144, 127] + har_c_lc = [10, 9, 8.5, 7, 6.5, 5, 4.5] + har_c_lc_angle = [230, 245, 160, 165, 155, 145, 135] + + harmonics = [1, 3, 5, 7, 9, 11, 13, 15] + + for i in range(1, 7): + pp.create_bus(net, 0.4, index=i) + + pp.create_ext_grid(net, bus=0, s_sc_max_mva=5, rx_max=3, r0x0_max=3, x0x_max=3) + + pp.create_transformer_from_parameters(net, hv_bus=0, lv_bus=1, sn_mva=0.4, vn_hv_kv=10, vn_lv_kv=0.4, vkr_percent=0, + vk_percent=4,pfe_kw=0, i0_percent=0, vector_group='YNyn',shift_degree=0, + vkr0_percent=0, vk0_percent=4, mag0_percent=100, mag0_rx=100, si0_hv_partial=0.9) + + pp.create_line_from_parameters(net, from_bus=1, to_bus=2, length_km=0.3, r_ohm_per_km=0.203, x_ohm_per_km=0.08, + c_nf_per_km=0, max_i_ka=275, r0_ohm_per_km=0.812, x0_ohm_per_km=0.24, c0_nf_per_km=0) + pp.create_line_from_parameters(net, from_bus=2, to_bus=3, length_km=0.1, r_ohm_per_km=0.606, x_ohm_per_km=0.083, + c_nf_per_km=0, max_i_ka=144, r0_ohm_per_km=2.424, x0_ohm_per_km=0.249, c0_nf_per_km=0) + pp.create_line_from_parameters(net, from_bus=2, to_bus=4, length_km=0.3, r_ohm_per_km=0.203, x_ohm_per_km=0.08, + c_nf_per_km=0, max_i_ka=275, r0_ohm_per_km=0.812, x0_ohm_per_km=0.24, c0_nf_per_km=0) + pp.create_line_from_parameters(net, from_bus=4, to_bus=5, length_km=0.1, r_ohm_per_km=0.606, x_ohm_per_km=0.083, + c_nf_per_km=0, max_i_ka=144, r0_ohm_per_km=2.424, x0_ohm_per_km=0.249, c0_nf_per_km=0) + pp.create_line_from_parameters(net, from_bus=4, to_bus=6, length_km=0.1, r_ohm_per_km=0.606, x_ohm_per_km=0.083, + c_nf_per_km=0, max_i_ka=144, r0_ohm_per_km=2.424, x0_ohm_per_km=0.249, c0_nf_per_km=0) + + pp.create_asymmetric_load(net, 3, p_a_mw=5/1000, p_b_mw=6/1000, p_c_mw=4/1000) + pp.create_asymmetric_load(net, 3, p_a_mw=0, p_b_mw=0, p_c_mw=0) + pp.create_asymmetric_load(net, 5, p_a_mw=6/1000, p_b_mw=7/1000, p_c_mw=8/1000) + pp.create_asymmetric_load(net, 5, p_a_mw=0, p_b_mw=0, p_c_mw=0) + pp.create_asymmetric_load(net, 6, p_a_mw=6/1000, p_b_mw=8/1000, p_c_mw=5/1000) + pp.create_asymmetric_load(net, 6, p_a_mw=0, p_b_mw=0, p_c_mw=0) + + # pp.runpp_3ph(net) + + # simple_plot(net, line_width=2.5, trafo_size=2.0, bus_size=1.5, ext_grid_size=1.5, line_color='black') + + pp_volt = [] + + for i in net.res_bus_3ph.index: + pp_volt.append(net.res_bus_3ph.vm_a_pu[i]*100) + pp_volt.append(net.res_bus_3ph.vm_b_pu[i]*100) + pp_volt.append(net.res_bus_3ph.vm_c_pu[i]*100) + + thd_a, thd_b, thd_c, har_vol_a, har_vol_b, har_vol_c = unbalanced_thd_voltage(net, harmonics, har_a,har_a_angle, + har_b, har_b_angle, har_c, har_c_angle, + har_a_lc, har_a_lc_angle, har_b_lc, + har_b_lc_angle, har_c_lc, har_c_lc_angle, + analysis_type="unbalanced") + + a, b, c, d = unbalanced_harmonic_current_voltage(net, harmonics, har_a, har_a_angle, har_b, har_b_angle, har_c, + har_c_angle, har_a_lc, har_a_lc_angle, har_b_lc, har_b_lc_angle, + har_c_lc, har_c_lc_angle, analysis_type="unbalanced") + + # y0=net._ppc0['internal']['Ybus'].todense() + # y1=net._ppc1['internal']['Ybus'].todense() + + +if __name__ == "__main__": + pytest.main([__file__, "-xs"])