nv_powerflow.py 4.35 KB
Newer Older
Markus Mirz's avatar
Markus Mirz committed
1
import numpy as np
2
3
4
import math
        
def Ymatrix_calc(branch, node):
Markus Mirz's avatar
Markus Mirz committed
5
    Ymatrix = np.zeros((node.num,node.num),dtype=np.complex)
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's avatar
Markus Mirz committed
23
24
25
    z = np.zeros(2*(node.num))
    h = np.zeros(2*(node.num))
    H = np.zeros((2*(node.num),2*(node.num)))
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's avatar
Markus Mirz committed
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)
45
            idx2 = idx1 + node.num
Markus Mirz's avatar
Markus Mirz committed
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])
50
51
        elif t == 'PV':
            z[m+1] = node.pwr_flow_values[2][i]
Markus Mirz's avatar
Markus Mirz committed
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)
55
            idx2 = idx1 + node.num
Markus Mirz's avatar
Markus Mirz committed
56
57
            H[m][idx1] = - np.real(Ymatrix[i][idx1])
            H[m][idx2] = np.imag(Ymatrix[i][idx1])
58
    
Markus Mirz's avatar
Markus Mirz committed
59
60
    epsilon = 10**(-10)
    diff = 5
Markus Mirz's avatar
Markus Mirz committed
61
    V = np.ones(node.num) + 1j* np.zeros(node.num)
62
63
    num_iter = 0
    
Markus Mirz's avatar
Markus Mirz committed
64
    State = np.ones(2*branch.num)
Markus Mirz's avatar
Markus Mirz committed
65
    State = np.concatenate((np.array([1,0]),State),axis=0)
66
    
Markus Mirz's avatar
Markus Mirz committed
67
    while diff > epsilon:
68
69
70
71
        for k in range(1,node.num+1):
            i = k-1
            m = 2*i
            i2 = i + node.num
Markus Mirz's avatar
Markus Mirz committed
72
            if node.type[i] == 'slack':
Markus Mirz's avatar
Markus Mirz committed
73
74
                h[m] = np.inner(H[m],State)
                h[m+1] = np.inner(H[m+1],State)
Markus Mirz's avatar
Markus Mirz committed
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's avatar
Markus Mirz committed
78
79
                h[m] = np.inner(H[m],State)
                h[m+1] = np.inner(H[m+1],State)
Markus Mirz's avatar
Markus Mirz committed
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's avatar
Markus Mirz committed
82
                h[m] = np.inner(H[m],State)
Markus Mirz's avatar
Markus Mirz committed
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]))
86
        
Markus Mirz's avatar
Markus Mirz committed
87
88
89
        r = np.subtract(z,h)
        Hinv = np.linalg.inv(H)
        Delta_State = np.inner(Hinv,r)
90
        State = State + Delta_State
Markus Mirz's avatar
Markus Mirz committed
91
        diff = np.amax(np.absolute(Delta_State))
92
        
Markus Mirz's avatar
Markus Mirz committed
93
        V = State[:node.num] + 1j * State[node.num:]
94
95
96
                
        num_iter = num_iter+1
        
Markus Mirz's avatar
Markus Mirz committed
97
    Irx = np.zeros((branch.num),dtype=np.complex)
98
99
100
    for idx in range(branch.num):
        fr = branch.start[idx]-1
        to = branch.end[idx]-1
Markus Mirz's avatar
Markus Mirz committed
101
102
103
        Irx[idx] = - (V[fr] - V[to])*Ymatrix[fr][to]
    Ir = np.real(Irx)
    Ix = np.imag(Irx)
104
    
Markus Mirz's avatar
Markus Mirz committed
105
106
107
    I = Ir + 1j*Ix
    Iinj_r = np.zeros(node.num)
    Iinj_x = np.zeros(node.num)
108
    for k in range(1,node.num+1):
Markus Mirz's avatar
Markus Mirz committed
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)
113
        
Markus Mirz's avatar
Markus Mirz committed
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's avatar
Markus Mirz committed
118
119
    S1 = np.multiply(V[branch.start-1], np.conj(I))
    S2 = - np.multiply(V[branch.end-1], np.conj(I))
120
121
122
    
    return V, I, Iinj, S1, S2, Sinj, num_iter