nv_powerflow.py 4.35 KB
 Markus Mirz committed Nov 07, 2018 1 ``````import numpy as np `````` Markus Mirz committed Nov 07, 2018 2 3 4 ``````import math def Ymatrix_calc(branch, node): `````` Markus Mirz committed Nov 07, 2018 5 `````` Ymatrix = np.zeros((node.num,node.num),dtype=np.complex) `````` Markus Mirz committed Nov 07, 2018 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 `````` Adjacencies = [[] for _ in range(node.num)] for index in range(branch.num): fr = branch.start[index] - 1 to = branch.end[index] - 1 Ymatrix[fr][to] = - branch.Y[index] Ymatrix[to][fr] = - branch.Y[index] Ymatrix[fr][fr] += branch.Y[index] Ymatrix[to][to] += branch.Y[index] Adjacencies[fr].append(to+1) Adjacencies[to].append(fr+1) return Ymatrix, Adjacencies def NV_power_flow(branch, node): """It performs Power Flow by using rectangular node voltage state variables.""" Ymatrix, Adj = Ymatrix_calc(branch,node) `````` Markus Mirz committed Nov 07, 2018 23 24 25 `````` z = np.zeros(2*(node.num)) h = np.zeros(2*(node.num)) H = np.zeros((2*(node.num),2*(node.num))) `````` Markus Mirz committed Nov 07, 2018 26 27 28 29 30 31 32 33 34 35 36 37 38 39 `````` for k in range(1,node.num+1): i = k-1 m = 2*i i2 = i + node.num t = node.type[i] if t == 'slack': V = node.pwr_flow_values[1][i] theta = node.pwr_flow_values[2][i] z[m] = V*math.cos(theta) z[m+1] = V*math.sin(theta) H[m][i] = 1 H[m+1][i2] = 1 elif t == 'PQ': `````` Markus Mirz committed Nov 07, 2018 40 41 42 43 44 `````` H[m][i] = - np.real(Ymatrix[i][i]) H[m][i2] = np.imag(Ymatrix[i][i]) H[m+1][i] = - np.imag(Ymatrix[i][i]) H[m+1][i2] = - np.real(Ymatrix[i][i]) idx1 = np.subtract(Adj[i],1) `````` Markus Mirz committed Nov 07, 2018 45 `````` idx2 = idx1 + node.num `````` Markus Mirz committed Nov 07, 2018 46 47 48 49 `````` H[m][idx1] = - np.real(Ymatrix[i][idx1]) H[m][idx2] = np.imag(Ymatrix[i][idx1]) H[m+1][idx1] = - np.imag(Ymatrix[i][idx1]) H[m+1][idx2] = - np.real(Ymatrix[i][idx1]) `````` Markus Mirz committed Nov 07, 2018 50 51 `````` elif t == 'PV': z[m+1] = node.pwr_flow_values[2][i] `````` Markus Mirz committed Nov 07, 2018 52 53 54 `````` H[m][i] = - np.real(Ymatrix[i][i]) H[m][i2] = np.imag(Ymatrix[i][i]) idx1 = np.subtract(Adj[i],1) `````` Markus Mirz committed Nov 07, 2018 55 `````` idx2 = idx1 + node.num `````` Markus Mirz committed Nov 07, 2018 56 57 `````` H[m][idx1] = - np.real(Ymatrix[i][idx1]) H[m][idx2] = np.imag(Ymatrix[i][idx1]) `````` Markus Mirz committed Nov 07, 2018 58 `````` `````` Markus Mirz committed Nov 07, 2018 59 60 `````` epsilon = 10**(-10) diff = 5 `````` Markus Mirz committed Nov 08, 2018 61 `````` V = np.ones(node.num) + 1j* np.zeros(node.num) `````` Markus Mirz committed Nov 07, 2018 62 63 `````` num_iter = 0 `````` Markus Mirz committed Nov 08, 2018 64 `````` State = np.ones(2*branch.num) `````` Markus Mirz committed Nov 07, 2018 65 `````` State = np.concatenate((np.array([1,0]),State),axis=0) `````` Markus Mirz committed Nov 07, 2018 66 `````` `````` Markus Mirz committed Nov 07, 2018 67 `````` while diff > epsilon: `````` Markus Mirz committed Nov 07, 2018 68 69 70 71 `````` for k in range(1,node.num+1): i = k-1 m = 2*i i2 = i + node.num `````` Markus Mirz committed Nov 08, 2018 72 `````` if node.type[i] == 'slack': `````` Markus Mirz committed Nov 07, 2018 73 74 `````` h[m] = np.inner(H[m],State) h[m+1] = np.inner(H[m+1],State) `````` Markus Mirz committed Nov 08, 2018 75 76 77 `````` elif node.type[i] == 'PQ': z[m] = (node.pwr_flow_values[1][i]*V[i].real + node.pwr_flow_values[2][i]*V[i].imag)/(np.abs(V[i])**2) z[m+1] = (node.pwr_flow_values[1][i]*V[i].imag - node.pwr_flow_values[2][i]*V[i].real)/(np.abs(V[i])**2) `````` Markus Mirz committed Nov 07, 2018 78 79 `````` h[m] = np.inner(H[m],State) h[m+1] = np.inner(H[m+1],State) `````` Markus Mirz committed Nov 08, 2018 80 81 `````` elif node.type[i] == 'PV': z[m] = (node.pwr_flow_values[1][i]*V[i].real + node.pwr_flow_values[2][i]*V[i].imag)/(np.abs(V[i])**2) `````` Markus Mirz committed Nov 07, 2018 82 `````` h[m] = np.inner(H[m],State) `````` Markus Mirz committed Nov 08, 2018 83 84 85 `````` 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])) `````` Markus Mirz committed Nov 07, 2018 86 `````` `````` Markus Mirz committed Nov 07, 2018 87 88 89 `````` r = np.subtract(z,h) Hinv = np.linalg.inv(H) Delta_State = np.inner(Hinv,r) `````` Markus Mirz committed Nov 07, 2018 90 `````` State = State + Delta_State `````` Markus Mirz committed Nov 07, 2018 91 `````` diff = np.amax(np.absolute(Delta_State)) `````` Markus Mirz committed Nov 07, 2018 92 `````` `````` Markus Mirz committed Nov 07, 2018 93 `````` V = State[:node.num] + 1j * State[node.num:] `````` Markus Mirz committed Nov 07, 2018 94 95 96 `````` num_iter = num_iter+1 `````` Markus Mirz committed Nov 07, 2018 97 `````` Irx = np.zeros((branch.num),dtype=np.complex) `````` Markus Mirz committed Nov 07, 2018 98 99 100 `````` for idx in range(branch.num): fr = branch.start[idx]-1 to = branch.end[idx]-1 `````` Markus Mirz committed Nov 07, 2018 101 102 103 `````` Irx[idx] = - (V[fr] - V[to])*Ymatrix[fr][to] Ir = np.real(Irx) Ix = np.imag(Irx) `````` Markus Mirz committed Nov 07, 2018 104 `````` `````` Markus Mirz committed Nov 07, 2018 105 106 107 `````` I = Ir + 1j*Ix Iinj_r = np.zeros(node.num) Iinj_x = np.zeros(node.num) `````` Markus Mirz committed Nov 07, 2018 108 `````` for k in range(1,node.num+1): `````` Markus Mirz committed Nov 07, 2018 109 110 111 112 `````` to = np.where(branch.end==k) fr = np.where(branch.start==k) Iinj_r[k-1] = np.sum(I[to[0]].real) - np.sum(I[fr[0]].real) Iinj_x[k-1] = np.sum(I[to[0]].imag) - np.sum(I[fr[0]].imag) `````` Markus Mirz committed Nov 07, 2018 113 `````` `````` Markus Mirz committed Nov 07, 2018 114 115 116 117 `````` Iinj = Iinj_r + 1j * Iinj_x Sinj_rx = np.multiply(V, np.conj(Iinj)) Sinj = np.real(Sinj_rx) + 1j * np.imag(Sinj_rx) `````` Markus Mirz committed Nov 08, 2018 118 119 `````` S1 = np.multiply(V[branch.start-1], np.conj(I)) S2 = - np.multiply(V[branch.end-1], np.conj(I)) `````` Markus Mirz committed Nov 07, 2018 120 121 122 `````` return V, I, Iinj, S1, S2, Sinj, num_iter ``````