nv_powerflow.py 4.99 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
import numpy
import math

class Real_to_all:
    def __init__(self,Ar, Ax):
        self.real = Ar
        self.imag = Ax
        self.complex = Ar + 1j*Ax
        self.mag = numpy.absolute(self.complex)
        self.phase = numpy.angle(self.complex)
        
def Ymatrix_calc(branch, node):
    Ymatrix = numpy.zeros((node.num,node.num),dtype=numpy.complex)
    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)
    
    z = numpy.zeros(2*(node.num))
    h = numpy.zeros(2*(node.num))
    H = numpy.zeros((2*(node.num),2*(node.num)))
    
    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':
            H[m][i] = - numpy.real(Ymatrix[i][i])
            H[m][i2] = numpy.imag(Ymatrix[i][i])
            H[m+1][i] = - numpy.imag(Ymatrix[i][i])
            H[m+1][i2] = - numpy.real(Ymatrix[i][i])
            idx1 = numpy.subtract(Adj[i],1)
            idx2 = idx1 + node.num
            H[m][idx1] = - numpy.real(Ymatrix[i][idx1])
            H[m][idx2] = numpy.imag(Ymatrix[i][idx1])
            H[m+1][idx1] = - numpy.imag(Ymatrix[i][idx1])
            H[m+1][idx2] = - numpy.real(Ymatrix[i][idx1])
        elif t == 'PV':
            z[m+1] = node.pwr_flow_values[2][i]
            H[m][i] = - numpy.real(Ymatrix[i][i])
            H[m][i2] = numpy.imag(Ymatrix[i][i])
            idx1 = numpy.subtract(Adj[i],1)
            idx2 = idx1 + node.num
            H[m][idx1] = - numpy.real(Ymatrix[i][idx1])
            H[m][idx2] = numpy.imag(Ymatrix[i][idx1])
    
    epsilon = 5
    Vr = numpy.ones(node.num)
    Vx = numpy.zeros(node.num)
    V = Real_to_all(Vr, Vx)
    num_iter = 0
    
    StateVr = numpy.ones(node.num)
    StateVx = numpy.zeros(node.num)
    State = numpy.concatenate((StateVr,StateVx),axis=0)
    
    while epsilon>10**(-10):
        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':
                h[m] = numpy.inner(H[m],State)
                h[m+1] = numpy.inner(H[m+1],State)
            elif t == 'PQ':
                z[m] = (node.pwr_flow_values[1][i]*V.real[i] + node.pwr_flow_values[2][i]*V.imag[i])/(V.mag[i]**2)
                z[m+1] = (node.pwr_flow_values[1][i]*V.imag[i] - node.pwr_flow_values[2][i]*V.real[i])/(V.mag[i]**2)
                h[m] = numpy.inner(H[m],State)
                h[m+1] = numpy.inner(H[m+1],State)
            elif t == 'PV':
                z[m] = (node.pwr_flow_values[1][i]*V.real[i] + node.pwr_flow_values[2][i]*V.imag[i])/(V.mag[i]**2)
                h[m] = numpy.inner(H[m],State)
                h[m+1] = V.mag[i]
                H[m+1][i] = numpy.cos(V.phase[i])
                H[m+1][i2] = numpy.sin(V.phase[i])
        
        r = numpy.subtract(z,h)
        Hinv = numpy.linalg.inv(H)
        Delta_State = numpy.inner(Hinv,r)
        State = State + Delta_State
        epsilon = numpy.amax(numpy.absolute(Delta_State))
        
        V.real = State[:node.num]
        V.imag = State[node.num:]
        V = Real_to_all(V.real, V.imag)
                
        num_iter = num_iter+1
        
    Irx = numpy.zeros((branch.num),dtype=numpy.complex)
    for idx in range(branch.num):
        fr = branch.start[idx]-1
        to = branch.end[idx]-1
        Irx[idx] = - (V.complex[fr] - V.complex[to])*Ymatrix[fr][to]
    Ir = numpy.real(Irx)
    Ix = numpy.imag(Irx)
    
    I = Real_to_all(Ir,Ix)
    Iinj_r = numpy.zeros(node.num)
    Iinj_x = numpy.zeros(node.num)
    for k in range(1,node.num+1):
        to = numpy.where(branch.end==k)
        fr = numpy.where(branch.start==k)
        Iinj_r[k-1] = numpy.sum(I.real[to[0]]) - numpy.sum(I.real[fr[0]])
        Iinj_x[k-1] = numpy.sum(I.imag[to[0]]) - numpy.sum(I.imag[fr[0]])
        
    Iinj = Real_to_all(Iinj_r,Iinj_x)
    Sinj_rx = numpy.multiply(V.complex,numpy.conj(Iinj.complex))
    Sinj = Real_to_all(numpy.real(Sinj_rx),numpy.imag(Sinj_rx))
    S1_rx = numpy.multiply(V.complex[branch.start-1],numpy.conj(I.complex))
    S2_rx = - numpy.multiply(V.complex[branch.end-1],numpy.conj(I.complex))
    S1 = Real_to_all(numpy.real(S1_rx),numpy.imag(S1_rx))
    S2 = Real_to_all(numpy.real(S2_rx),numpy.imag(S2_rx))
    
    return V, I, Iinj, S1, S2, Sinj, num_iter