nv_powerflow.py 4.23 KB
Newer Older
 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 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 `````` 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: H[m][i] = - np.real(system.Ymatrix[i][i]) H[m][i2] = np.imag(system.Ymatrix[i][i]) H[m + 1][i] = - np.imag(system.Ymatrix[i][i]) H[m + 1][i2] = - np.real(system.Ymatrix[i][i]) idx1 = np.subtract(system.Adjacencies[i], 1) idx2 = idx1 + nodes_num H[m][idx1] = - np.real(system.Ymatrix[i][idx1]) H[m][idx2] = np.imag(system.Ymatrix[i][idx1]) H[m + 1][idx1] = - np.imag(system.Ymatrix[i][idx1]) H[m + 1][idx2] = - np.real(system.Ymatrix[i][idx1]) elif node_type is BusType.PV: z[m + 1] = np.real(node.power) H[m][i] = - np.real(system.Ymatrix[i][i]) H[m][i2] = np.imag(system.Ymatrix[i][i]) idx1 = np.subtract(system.Adjacencies[i], 1) idx2 = idx1 + nodes_num H[m][idx1] = - np.real(system.Ymatrix[i][idx1]) H[m][idx2] = np.imag(system.Ymatrix[i][idx1]) `````` Jan Dinkelbach committed Oct 25, 2019 61 62 `````` epsilon = 10 ** (-10) `````` martin.moraga committed Dec 11, 2019 63 `````` #epsilon = 0.01 `````` Jan Dinkelbach committed Oct 25, 2019 64 65 66 67 68 69 70 `````` 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 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 `````` 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]) + np.imag(node.power_pu) * np.imag(V[i]))(np.abs(V[i]) ** 2) 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 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 `````` 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``````