Commit 3135af6f authored by Markus Mirz's avatar Markus Mirz
Browse files

clean up powerflow module

parent 8b890869
......@@ -2,7 +2,7 @@ import numpy
import math
import matplotlib as plt
import Network_data
import Power_flow
import powerflow
import Meas_data
import StateEstimator
......@@ -50,7 +50,7 @@ slackV = 1.02
Base = PerUnit(S,V)
branch, node = Network_data.Network_95_nodes(Base, slackV)
Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue, num_iter = Power_flow.BC_power_flow(branch, node)
Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue, num_iter = powerflow.BC_power_flow(branch, node)
""" Write here the indexes of the nodes/branches where the measurements are"""
V_idx = numpy.array([1,11,55])
......
def BC_power_flow(branch, node):
"""It performs Power Flow by using branch current state variables."""
import numpy
import math
class Znod:
def __init__(self,branch,node):
self.Z = numpy.zeros((node.num-1,branch.num),dtype=numpy.complex_)
for k in range(1,node.num):
i = k-1
aa = branch.start[i]-2
if aa == -1:
self.Z[i] = self.Z[i]
else:
self.Z[i] = self.Z[aa]
self.Z[i][i] = branch.Z[i]
zzz = numpy.zeros((1,branch.num))
self.Z = numpy.concatenate((zzz,self.Z), axis=0)
self.R = self.Z.real
self.X = self.Z.imag
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)
znod = Znod(branch,node)
z = numpy.zeros(2*(branch.num+1))
h = numpy.zeros(2*(branch.num+1))
H = numpy.zeros((2*(branch.num+1),2*(branch.num+1)))
for k in range(1,node.num+1):
i = k-1
m = 2*i
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][0] = 1;
H[m][2:] = numpy.concatenate((-znod.R[i],znod.X[i]),axis=0)
H[m+1][1] = 1
H[m+1][2:] = numpy.concatenate((-znod.X[i],-znod.R[i]),axis=0)
elif t == 'PQ':
to = numpy.where(branch.end==k)
fr = numpy.where(branch.start==k)
to1 = to[0]+2
to2 = to1+branch.num
fr1 = fr[0]+2
fr2 = fr1+branch.num
H[m][to1] = 1
H[m+1][to2] = 1
H[m][fr1] = -1
H[m+1][fr2] = -1
elif t == 'PV':
z[m+1] = node.pwr_flow_values[2][i]
to = numpy.where(branch.end==k)
fr = numpy.where(branch.start==k)
to1 = to[0]+2
fr1 = fr[0]+2
H[m][to1] = 1
H[m][fr1] = -1
epsilon = 5
Vr = numpy.ones(node.num)
Vx = numpy.zeros(node.num)
V = Real_to_all(Vr, Vx)
num_iter = 0
State = numpy.zeros(2*branch.num)
State = numpy.concatenate((numpy.array([1,0]),State),axis=0)
while epsilon>10**(-10):
for k in range(1,node.num+1):
i = k-1
m = 2*i
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][0] = numpy.cos(V.theta[i])
H[m+1][1] = numpy.sin(V.theta[i])
idx = numpy.where(znod.Z[i] != 0)
idx1 = idx[0]+2
idx2 = idx1[0]+branch.num
H[m+1][idx1] = - znod.R[i][idx]*numpy.cos(V.theta[i]) - znod.X[i][idx]*numpy.sin(V.theta[i])
H[m+1][idx2] = - znod.R[i][idx]*numpy.sin(V.theta[i]) - znod.X[i][idx]*numpy.cos(V.theta[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[0] = State[0]
V.imag[0] = State[1]
Ir = State[2:branch.num+2]
Ix = State[branch.num+2:]
Irx = Ir + 1j*Ix
Vcomplex = numpy.inner((V.real[0] + 1j*V.imag[0]),numpy.ones(node.num))
DeltaV = numpy.inner(znod.Z,Irx)
V.complex = Vcomplex[0] - DeltaV
V.real = numpy.real(V.complex)
V.imag = numpy.imag(V.complex)
V = Real_to_all(V.real, V.imag)
num_iter = num_iter+1
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
\ No newline at end of file
import numpy as np
import math
class Znod:
def __init__(self, branches, nodes):
# Complex impedance matrix
self.Z = np.zeros( (nodes.num-1,branches.num), dtype=np.complex_ )
# Add impedance of each branch
for k in range(1, nodes.num):
i = k-1
aa = branches.start[i]-2
if aa == -1:
self.Z[i] = self.Z[i]
else:
self.Z[i] = self.Z[aa]
self.Z[i][i] = branches.Z[i]
# Add one row of zeros
zzz = np.zeros((1, branches.num))
self.Z = np.concatenate( (zzz,self.Z), axis=0 )
# Resistance matrix
self.R = self.Z.real
# Reactance matrix
self.X = self.Z.imag
def BC_power_flow(branch, node):
"""It performs Power Flow by using branch current state variables."""
znod = Znod(branch,node)
# real matrices storing complex numbers
# real in pos n and imag in pos n+1
z = np.zeros(2*(branch.num+1))
h = np.zeros(2*(branch.num+1))
H = np.zeros((2*(branch.num+1),2*(branch.num+1)))
for k in range(1, node.num+1):
# i is matrix index
# the size of z,h,H is twice the number of branches
# but here we are using the node number to work on the
# matrix. Is this not dangerous?
i = k-1
m = 2*i
if node.type[i] == 'slack':
V = node.pwr_flow_values[1][i]
theta = node.pwr_flow_values[2][i]
# real node voltage
z[m] = V*math.cos(theta)
# imag node voltage
z[m+1] = V*math.sin(theta)
H[m][0] = 1;
H[m][2:] = np.concatenate((-znod.R[i],znod.X[i]), axis=0)
H[m+1][1] = 1
H[m+1][2:] = np.concatenate((-znod.X[i],-znod.R[i]), axis=0)
elif node.type[i] == 'PQ':
to = np.where(branch.end==k)
fr = np.where(branch.start==k)
to1 = to[0]+2
to2 = to1+branch.num
fr1 = fr[0]+2
fr2 = fr1+branch.num
H[m][to1] = 1
H[m+1][to2] = 1
H[m][fr1] = -1
H[m+1][fr2] = -1
elif node.type[i] == 'PV':
z[m+1] = node.pwr_flow_values[2][i]
to = np.where(branch.end==k)
fr = np.where(branch.start==k)
to1 = to[0]+2
fr1 = fr[0]+2
H[m][to1] = 1
H[m][fr1] = -1
print(z)
print(h)
print(H)
epsilon = 5
V = np.ones(node.num) + 1j* np.zeros(node.num)
num_iter = 0
State = np.zeros(2*branch.num)
State = np.concatenate( (np.array([1,0]),State), axis=0 )
while epsilon > 10**(-10):
for k in range(1,node.num+1):
i = k-1
m = 2*i
t = node.type[i]
if t == 'slack':
h[m] = np.inner(H[m],State)
h[m+1] = np.inner(H[m+1],State)
elif t == '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)
h[m] = np.inner(H[m],State)
h[m+1] = np.inner(H[m+1],State)
elif t == '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)
h[m] = np.inner(H[m],State)
h[m+1] = np.abs(V[i])
H[m+1][0] = np.cos(np.angle(V[i]))
H[m+1][1] = np.sin(np.angle(V[i]))
idx = np.where(znod.Z[i] != 0)
idx1 = idx[0]+2
idx2 = idx1[0]+branch.num
H[m+1][idx1] = - znod.R[i][idx]*np.cos(np.angle(V[i])) - znod.X[i][idx]*np.sin(np.angle(V[i]))
H[m+1][idx2] = - znod.R[i][idx]*np.sin(np.angle(V[i])) - znod.X[i][idx]*np.cos(np.angle(V[i]))
r = np.subtract(z,h)
Hinv = np.linalg.inv(H)
Delta_State = np.inner(Hinv,r)
State = State + Delta_State
epsilon = np.amax(np.absolute(Delta_State))
V[0] = State[0] + 1j * State[1]
Ir = State[2:branch.num+2]
Ix = State[branch.num+2:]
Irx = Ir + 1j*Ix
Vcomplex = np.inner((V[0].real + 1j*V[0].imag), np.ones(node.num))
DeltaV = np.inner(znod.Z,Irx)
V = Vcomplex[0] - DeltaV
num_iter = num_iter+1
I = Ir + 1j*Ix
Iinj_r = np.zeros(node.num)
Iinj_x = np.zeros(node.num)
for k in range(1,node.num+1):
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)
Iinj = Iinj_r + 1j * Iinj_x
Sinj_rx = np.multiply(V, np.conj(Iinj))
Sinj = np.real(Sinj_rx) + 1j * np.imag(Sinj_rx)
S1_rx = np.multiply(V[branch.start-1], np.conj(I))
S2_rx = np.multiply(V[branch.end-1], np.conj(I))
S1 = np.real(S1_rx) + 1j * np.imag(S1_rx)
S2 = np.real(S2_rx) + 1j * np.imag(S2_rx)
return V, I, Iinj, S1, S2, Sinj, num_iter
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment