Commit e1217015 authored by Markus Mirz's avatar Markus Mirz
Browse files

Merge branch 'restructure' into 'master'

Restructure

See merge request acs/core/simulation/state-estimation!2
parents 25f6ac20 5cb610b2
......@@ -4,7 +4,8 @@ import matplotlib as plt
import Network_data
import Power_flow
import Meas_data
import StateEstimator
import NvStateEstimator
#import BcStateEstimator
class PerUnit:
def __init__(self, S, V):
......@@ -36,7 +37,6 @@ class Zdata_init:
nmeas = meas.V.num + meas.I.num + 2*meas.Sinj.num + 2*meas.S1.num + 2*meas.S2.num + meas.Ipmu_mag.num + meas.Ipmu_phase.num + meas.Vpmu_mag.num + meas.Vpmu_phase.num
self.mtype = numpy.zeros(nmeas)
self.mval = numpy.zeros(nmeas)
self.mmeas = numpy.zeros(nmeas)
self.mbranch = numpy.zeros(nmeas)
self.mfrom = numpy.zeros(nmeas)
self.mto = numpy.zeros(nmeas)
......@@ -50,20 +50,20 @@ 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 = Power_flow.NV_power_flow(branch, node)
""" Write here the indexes of the nodes/branches where the measurements are"""
V_idx = numpy.array([1,11,55])
I_idx = numpy.array([])
I_idx = numpy.array([13,36])
Sinj_idx = numpy.linspace(2,node.num,node.num-1)
S1_idx = numpy.array([1,11,28,55,59])
S2_idx = numpy.array([10,54])
Ipmu_idx = numpy.array([])
Vpmu_idx = numpy.array([])
Ipmu_idx = numpy.array([1,11,55])
Vpmu_idx = numpy.array([13,36])
""" Write here the percent uncertainties of the measurements"""
V_unc = 1
I_unc = 0
I_unc = 2
Sinj_unc = 5
S_unc = 2
Pmu_mag_unc = 0.7
......@@ -93,13 +93,15 @@ MC_trials = 1
#Sinjest_matrix = numpy.zeros((MC_trials,node.num),dtype=numpy.complex_)
iter_counter = numpy.zeros(MC_trials)
Ymatr, Adj = Power_flow.Ymatrix_calc(branch,node)
for mciter in range(0,MC_trials):
if mciter%10 == 0:
print(mciter)
zdatameas = Meas_data.Zdatameas_creation(zdata, zdatameas)
Vest, Iest, Iinjest, S1est, S2est, Sinjest = StateEstimator.BcDsseCall(branch, node, zdatameas)
Vest, Iest, Iinjest, S1est, S2est, Sinjest = NvStateEstimator.DsseCall(branch, node, zdatameas, Ymatr, Adj)
# Iest_matrix[mciter] = Iest.complex
# Vest_matrix[mciter] = Vest.complex
......
......@@ -12,6 +12,7 @@ def Network_95_nodes(Base, slackV):
# self.R = numpy.divide(R,Base.Z)
# self.X = numpy.divide(X,Base.Z)
self.Z = self.R + 1j*self.X
self.Y = 1/self.Z
class Node:
def __init__(self, ntype, slackV, P, Q):
......@@ -57,20 +58,20 @@ def Network_95_nodes(Base, slackV):
nPQ = 'PQ'
P = [0, 0.82092, 0, 0, 0.01461, 0, 0.03046, 0, 0.02805, 0, 0.00487, 0, 0.00487, 0, 0.01569, 0, 0.01316, 0, 0,
P = numpy.array([0, 0.82092, 0, 0, 0.01461, 0, 0.03046, 0, 0.02805, 0, 0.00487, 0, 0.00487, 0, 0.01569, 0, 0.01316, 0, 0,
0.05754, 0.05725, 0.02777, 0.0292, 0, 0.01464, 0.6776, -0.675, 0, 0.02632, 0, 0.01831, 0, 0.0483, 0.03407,
0.02777, 0, 0.05605, 0.11734, 0.00733, 0.01098, 0, 0, 0.03167, 0.05264, 0, 0, 0, 0.01461, 0, 0.01562, 0,
0.03167, 0.04493, 0, 0, 0.00779, 0.00876, 0.00733, 0.02777, 0, 0.01461, 0.01461, 0.01783, 0, 0.02681, 0.0292,
0.02734, 0, 0.03407, 0.01363, 0, 0, 0, 0.01266, 0.05078, 0, 0, 0, 0, 0.05359, 0, 0.025, 0, 0.08578, 0, 0.09698,
0.04563, 0, 0.13277, 0, 0.0561, 0, 0.04806, -0.675]
P = numpy.multiply(P,(10**6/Base.S))
Q = [0, 0.14346, 0, 0, 0.00208, 0, 0.00434, 0, 0.004, 0, 0.00069, 0, 0.00069, 0, 0.00516, 0, 0.00224, 0, 0, 0.00969,
0.04563, 0, 0.13277, 0, 0.0561, 0, 0.04806, -0.675])
P = P*(10**6/Base.S)
Q = numpy.array([0, 0.14346, 0, 0, 0.00208, 0, 0.00434, 0, 0.004, 0, 0.00069, 0, 0.00069, 0, 0.00516, 0, 0.00224, 0, 0, 0.00969,
0.00817, 0.00432, 0.00416, 0, 0.00209, 0.22267, 0.2219, 0, 0.00449, 0, 0.00261, 0, 0.00762, 0.00486, 0.00432,
0, 0.00799, 0.02783, 0.00104, 0.00156, 0, 0, 0.00451, 0.00899, 0, 0, 0, 0.00208, 0, 0.00259, 0, 0.00451,
0.00863, 0, 0, 0.00111, 0.00124, 0.00104, 0.00432, 0, 0.00208, 0.00208, 0.00328, 0, 0.00381, 0.00416, 0.005,
0, 0.00486, 0.00194, 0, 0, 0, 0.00181, 0.00982, 0, 0, 0, 0, 0.00764, 0, 0.005, 0, 0.01897, 0, 0.01456, 0.01,
0, 0.02485, 0, 0.00799, 0, 0.00796, 0.2219]
Q = numpy.multiply(Q,(10**6/Base.S))
0, 0.02485, 0, 0.00799, 0, 0.00796, 0.2219])
Q = Q*(10**6/Base.S)
ntype = ['slack']
for j in range(1, len(P)+1):
......
import numpy
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 __init__(self,Ar, Ax):
""" Converts real and imaginary parts to all the other formats"""
self.real = Ar
self.imag = Ax
self.complex = Ar + 1j*Ax
self.mag = numpy.absolute(self.complex)
self.phase = numpy.angle(self.complex)
class Znod:
def __init__(self,branch,node):
""" Calculates a matrix of impedances for calculating the voltage in all
the nodes of the grid, starting from the slack bus voltage and knowing
the branch currents.
In each row it gives the set of branch impedances to calculate the voltage
of the corresponding node through the voltage drops starting from the
first node of the grid """
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
def BcDsseCall(branch, node, zdata):
def DsseCall(branch, node, zdata):
""" It identifies the type of measurements present in the measurement set and
calls the appropriate estimator for dealing with them."""
......@@ -26,12 +51,13 @@ def BcDsseCall(branch, node, zdata):
est_code = trad_code + PMU_code
if est_code == 1:
Vest, Iest = BcDsseTrad(branch, node, zdata)
Vest, Iest = DsseTrad(branch, node, zdata)
elif est_code == 2:
Vest, Iest = BcDssePmu(branch, node, zdata)
Vest, Iest = DssePmu(branch, node, zdata)
else:
Vest, Iest = BcDsseMixed(branch, node, zdata)
Vest, Iest = DsseMixed(branch, node, zdata)
""" From here on, all the other quantities of the grid are calculated. """
Iinj_r = numpy.zeros(node.num)
Iinj_x = numpy.zeros(node.num)
for k in range(1,node.num+1):
......@@ -44,34 +70,17 @@ def BcDsseCall(branch, node, zdata):
Sinj_rx = numpy.multiply(Vest.complex,numpy.conj(Iinjest.complex))
Sinjest = Real_to_all(numpy.real(Sinj_rx),numpy.imag(Sinj_rx))
S1_rx = numpy.multiply(Vest.complex[branch.start-1],numpy.conj(Iest.complex))
S2_rx = numpy.multiply(Vest.complex[branch.end-1],numpy.conj(Iest.complex))
S2_rx = -numpy.multiply(Vest.complex[branch.end-1],numpy.conj(Iest.complex))
S1est = Real_to_all(numpy.real(S1_rx),numpy.imag(S1_rx))
S2est = Real_to_all(numpy.real(S2_rx),numpy.imag(S2_rx))
return Vest, Iest, Iinjest, S1est, S2est, Sinjest
def BcDsseTrad(branch, node, zdata):
def DsseTrad(branch, node, zdata):
""" It performs state estimation using branch current state variables and it is customized
to work without PMU measurements"""
to work without PMU measurements. """
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
vidx = numpy.where(zdata.mtype==1)
pidx = numpy.where(zdata.mtype==2)
qidx = numpy.where(zdata.mtype==3)
......@@ -88,7 +97,7 @@ def BcDsseTrad(branch, node, zdata):
busvi = zdata.mfrom[vidx]
buspi = zdata.mfrom[pidx]
busqi = zdata.mfrom[qidx]
busqi = zdata.mfrom[qidx]
buspf = zdata.mfrom[pfidx]
busqf = zdata.mfrom[qfidx]
branchpf = zdata.mbranch[pfidx]
......@@ -101,8 +110,8 @@ def BcDsseTrad(branch, node, zdata):
Pinj = z[pidx]
Qinj = z[qidx]
Pbr = z[pfidx]*numpy.sign(branchpf)
Qbr = z[qfidx]*numpy.sign(branchqf)
Pbr = z[pfidx]
Qbr = z[qfidx]
idx = numpy.where(zdata.mstddev<10**(-6))
zdata.mstddev[idx] = 10**(-6)
......@@ -111,6 +120,8 @@ def BcDsseTrad(branch, node, zdata):
znod = Znod(branch,node)
""" Jacobian for Power Injection Measurements (converted to equivalent
rectangualar current measurements) """
H2 = numpy.zeros((npi,2*branch.num+1))
H3 = numpy.zeros((nqi,2*branch.num+1))
for i in range(0,npi):
......@@ -126,6 +137,8 @@ def BcDsseTrad(branch, node, zdata):
H3[i][to2] = 1
H3[i][fr2] = -1
""" Jacobian for branch Power Measurements (converted to equivalent
rectangualar current measurements)"""
H4 = numpy.zeros((npf,2*branch.num+1))
H5 = numpy.zeros((nqf,2*branch.num+1))
for i in range(0,npf):
......@@ -146,6 +159,7 @@ def BcDsseTrad(branch, node, zdata):
State = numpy.concatenate((numpy.array([1,]),State),axis=0)
while epsilon>10**(-6):
""" Computation of equivalent current measurements in place of the power measurements """
Irinj = (Pinj*V.real[buspi-1] + Qinj*V.imag[busqi-1])/(V.mag[buspi-1]**2)
Ixinj = (Pinj*V.imag[buspi-1] - Qinj*V.real[busqi-1])/(V.mag[buspi-1]**2)
z[pidx] = Irinj
......@@ -156,6 +170,7 @@ def BcDsseTrad(branch, node, zdata):
z[pfidx] = Irbr
z[qfidx] = Ixbr
""" Voltage Magnitude Measurements """
h1 = V.mag[busvi-1]
H1 = numpy.zeros((nvi,2*branch.num+1))
for i in range(0,nvi):
......@@ -168,12 +183,15 @@ def BcDsseTrad(branch, node, zdata):
H1[i][idx[0]+1] = -znod.R[m][idx]*numpy.cos(V.phase[m]) - znod.X[m][idx]*numpy.sin(V.phase[m])
H1[i][idx[0]+1+branch.num] = -znod.R[m][idx]*numpy.sin(V.phase[m]) + znod.X[m][idx]*numpy.cos(V.phase[m])
""" Power Injection Measurements """
h2 = numpy.inner(H2,State)
h3 = numpy.inner(H3,State)
""" Power Flow Measurements """
h4 = numpy.inner(H4,State)
h5 = numpy.inner(H5,State)
""" Current Magnitude Measurements """
h6 = I.mag[branchii]
phi = I.phase[branchii]
H6 = numpy.zeros((nii,2*branch.num+1))
......@@ -182,6 +200,7 @@ def BcDsseTrad(branch, node, zdata):
H6[i][idx] = numpy.cos(phi[i])
H6[i][idx+branch.num] = numpy.sin(phi[i])
""" WLS computation """
H = numpy.concatenate((H1,H2,H3,H4,H5,H6),axis=0)
y = numpy.concatenate((h1,h2,h3,h4,h5,h6),axis=0)
res = numpy.subtract(z,y)
......@@ -200,6 +219,7 @@ def BcDsseTrad(branch, node, zdata):
Ix = State[branch.num+1:]
Irx = Ir + 1j*Ix
""" Forward sweep to calculate node voltage along the whole grid """
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
......@@ -214,7 +234,7 @@ def BcDsseTrad(branch, node, zdata):
def BcDssePmu(branch, node, zdata):
def DssePmu(branch, node, zdata):
""" It performs state estimation using branch current state variables and it is customized
to work using only PMU measurements."""
......@@ -238,10 +258,9 @@ def BcDssePmu(branch, node, zdata):
def BcDsseMixed(branch, node, zdata):
def DsseMixed(branch, node, zdata):
""" It performs state estimation using branch current state variables and it is built
to work in scenarios where both conventional and PMU measurements are simultaneously present."""
import numpy
class Znod:
def __init__(self,branch,node):
......
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 BC_power_flow(branch, node):
"""It performs Power Flow by using branch current state variables."""
import numpy
import math
"""It performs Power Flow by using rectangular branch current state variables."""
class Znod:
def __init__(self,branch,node):
self.Z = numpy.zeros((node.num-1,branch.num),dtype=numpy.complex_)
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
......@@ -19,15 +43,6 @@ def BC_power_flow(branch, node):
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))
......@@ -92,13 +107,13 @@ def BC_power_flow(branch, node):
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])
H[m+1][0] = numpy.cos(V.phase[i])
H[m+1][1] = numpy.sin(V.phase[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])
H[m+1][idx1] = - znod.R[i][idx]*numpy.cos(V.phase[i]) - znod.X[i][idx]*numpy.sin(V.phase[i])
H[m+1][idx2] = - znod.R[i][idx]*numpy.sin(V.phase[i]) - znod.X[i][idx]*numpy.cos(V.phase[i])
r = numpy.subtract(z,h)
Hinv = numpy.linalg.inv(H)
......@@ -134,8 +149,122 @@ def BC_power_flow(branch, node):
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))
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
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
\ No newline at end of file
return V, I, Iinj, S1, S2, Sinj, num_iter
\ No newline at end of file
This diff is collapsed.
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][