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..d26d4f01c0 --- /dev/null +++ b/pandapower/harmonics/balanced.py @@ -0,0 +1,159 @@ +import math +import cmath +import numpy as np +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: 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) + + 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: 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 = [] + 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: + 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..206f3b5d09 --- /dev/null +++ b/pandapower/harmonics/harmonic_impedance_creator.py @@ -0,0 +1,196 @@ +# Script for creating impedance matrix for higher order frequencies + +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: pandapowerNet, harmon_order: list[int], analysis_type: str): + 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 = makeYbus(ppci_0["baseMVA"], ppci_0["bus"], ppci_0["branch"]) + + ppc_1, ppci_1 = _pd2ppc(net, 1) + + 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 = 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..9ddef24249 --- /dev/null +++ b/pandapower/harmonics/unbalanced.py @@ -0,0 +1,322 @@ +import math +import cmath +import numpy as np +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 + + +# 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(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 = [] + 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]) / (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 + 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) / 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(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 / net.res_bus_3ph.vm_a_pu[i + 1]) + harmonics_voltage_b[i + 1, j] = harmonics_voltage[3 * i + 1, j] * 100 * ( + 1 / net.res_bus_3ph.vm_b_pu[i + 1]) + harmonics_voltage_c[i + 1, j] = harmonics_voltage[3 * i + 2, j] * 100 * ( + 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 / 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(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(net, harmonics, analysis_type) + + s_base = net.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: + + runpp_3ph(net) + + current = [] + current_res = [] + current_pv = [] + current_0 = [] + + # p.u. values + 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 != 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 net.asymmetric_load.bus[j] == i: + connected = 1 # an asymmetric load is connected to the observed node + + ##Residential + 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(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(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(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(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(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(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(net.res_bus_3ph.vm_c_pu[i], + net.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(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(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(net.res_bus_3ph.vm_c_pu[i], + net.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(net.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/pd2ppc_zero.py b/pandapower/pd2ppc_zero.py index 3fd70cb4fd..0ef02ed661 100644 --- a/pandapower/pd2ppc_zero.py +++ b/pandapower/pd2ppc_zero.py @@ -254,10 +254,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/__init__.py b/pandapower/test/harmonics/__init__.py new file mode 100644 index 0000000000..e69de29bb2 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"]) 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 0000000000..6ced212d2c Binary files /dev/null and b/pandapower/test/test_files/harmonics/CIGRE_LV.xlsx differ diff --git a/pandapower/test/test_files/harmonics/Test_1f_Loads.xlsx b/pandapower/test/test_files/harmonics/Test_1f_Loads.xlsx new file mode 100644 index 0000000000..a7fd0aef81 Binary files /dev/null and b/pandapower/test/test_files/harmonics/Test_1f_Loads.xlsx differ 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 0000000000..2cedd8b273 Binary files /dev/null and b/pandapower/test/test_files/harmonics/Test_1f_lines.xlsx differ