nv_powerflow.py 3.65 KB
 Markus Mirz committed Nov 07, 2018 1 ``````import numpy as np `````` Jan Dinkelbach committed May 08, 2019 2 3 ``````from .network import BusType from .results import Results `````` Markus Mirz committed Nov 07, 2018 4 `````` `````` Jan Dinkelbach committed Oct 25, 2019 5 `````` `````` Jan Dinkelbach committed May 08, 2019 6 ``````def solve(system): `````` Jan Dinkelbach committed Oct 25, 2019 7 `````` """It performs powerflow by using rectangular node voltage state variables and considering the current mismatch function. `````` Jan Dinkelbach committed Oct 21, 2019 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 `````` Solve the non-linear powerflow problem stated by r = z-h(state) = 0 following the Newton-Raphson approach delta_state = H^-1 * r new_state = old_state + delta_state r: residual function (current mismatch) z: expected currents state: rectangular voltages (i.e. [V0_re, V1_re, ..., VN_re, V0_im, V1_im, ... , VN_im]) h: currents calculated from state H: Jacobian matrix V: same as state but with complex numbers (i.e. [V0_re+j*V0_im, V1_re+j*V1_im, ...]) """ `````` Jan Dinkelbach committed Oct 25, 2019 25 `````` `````` martin.moraga committed Dec 09, 2019 26 `````` nodes_num = system.get_nodes_num() `````` Jan Dinkelbach committed Oct 25, 2019 27 28 29 30 `````` z = np.zeros(2 * nodes_num) h = np.zeros(2 * nodes_num) H = np.zeros((2 * nodes_num, 2 * nodes_num)) `````` martin.moraga committed Dec 09, 2019 31 32 33 34 35 36 37 38 39 40 41 42 `````` for node in system.nodes: if node.ideal_connected_with =='': i = node.index m = 2 * i i2 = i + nodes_num node_type = node.type if node_type == BusType.SLACK: z[m] = np.real(node.voltage_pu) z[m + 1] = np.imag(node.voltage_pu) H[m][i] = 1 H[m + 1][i2] = 1 elif node_type is BusType.PQ: `````` admin-mpa committed Feb 11, 2020 43 44 45 46 `````` H[m][:node.num] = numpy.real(Ymatrix[i]) H[m][node.num:] = - numpy.imag(Ymatrix[i]) H[m+1][:node.num] = numpy.imag(Ymatrix[i]) H[m+1][node.num:] = numpy.real(Ymatrix[i]) `````` martin.moraga committed Dec 09, 2019 47 48 `````` elif node_type is BusType.PV: z[m + 1] = np.real(node.power) `````` admin-mpa committed Feb 11, 2020 49 50 `````` H[m][:node.num] = numpy.real(Ymatrix[i]) H[m][node.num:] = - numpy.imag(Ymatrix[i]) `````` Jan Dinkelbach committed Oct 25, 2019 51 52 `````` epsilon = 10 ** (-10) `````` martin.moraga committed Dec 11, 2019 53 `````` #epsilon = 0.01 `````` Jan Dinkelbach committed Oct 25, 2019 54 55 56 57 58 59 60 `````` diff = 5 V = np.ones(nodes_num) + 1j * np.zeros(nodes_num) num_iter = 0 state = np.concatenate((np.ones(nodes_num), np.zeros(nodes_num)), axis=0) while diff > epsilon: `````` martin.moraga committed Dec 09, 2019 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 `````` for node in system.nodes: if node.ideal_connected_with =='': i = node.index m = 2 * i i2 = i + nodes_num node_type = node.type if node_type is BusType.SLACK: h[m] = np.inner(H[m], state) h[m + 1] = np.inner(H[m + 1], state) elif node_type is BusType.PQ: z[m] = (np.real(node.power_pu) * np.real(V[i]) + np.imag(node.power_pu) * np.imag(V[i])) / (np.abs(V[i]) ** 2) z[m + 1] = (np.real(node.power_pu) * np.imag(V[i]) - np.imag(node.power_pu) * np.real(V[i])) / (np.abs(V[i]) ** 2) h[m] = np.inner(H[m], state) h[m + 1] = np.inner(H[m + 1], state) elif node_type is BusType.PV: z[m] = (np.real(node.power_pu) * np.real(V[i]) + `````` admin-mpa committed Feb 10, 2020 79 `````` np.imag(node.power_pu) * np.imag(V[i])) / (np.abs(V[i]) ** 2) `````` martin.moraga committed Dec 09, 2019 80 81 82 83 `````` h[m] = np.inner(H[m], state) h[m + 1] = np.abs(V[i]) H[m + 1][i] = np.cos(np.angle(V[i])) H[m + 1][i2] = np.sin(np.angle(V[i])) `````` Jan Dinkelbach committed Oct 25, 2019 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 `````` r = np.subtract(z, h) Hinv = np.linalg.inv(H) delta_state = np.inner(Hinv, r) state = state + delta_state diff = np.amax(np.absolute(delta_state)) V = state[:nodes_num] + 1j * state[nodes_num:] num_iter = num_iter + 1 # calculate all the other quantities of the grid powerflow_results = Results(system) powerflow_results.load_voltages(V) powerflow_results.calculate_all() return powerflow_results, num_iter``````