Commit 8c8f8f1b authored by Markus Mirz's avatar Markus Mirz
Browse files

initial commit of single phase estimator

parent 23ffd6fa
import numpy
import math
import matplotlib as plt
import Network_data
import Power_flow
import Meas_data
import StateEstimator
class PerUnit:
def __init__(self, S, V):
self.S = S
self.V = V
self.I = S/V
self.Z = S/(V**2)
class Measurements:
def __init__(self, index, unc):
self.index = index.astype(int)
self.unc = unc
self.num = len(index)
class Measurement_set:
def __init__(self, V, I, Sinj, S1, S2, Vpmu_mag, Vpmu_phase, Ipmu_mag, Ipmu_phase):
self.V = V
self.Sinj = Sinj
self.S1 = S1
self.S2 = S2
self.I = I
self.Vpmu_mag = Vpmu_mag
self.Vpmu_phase = Vpmu_phase
self.Ipmu_mag = Ipmu_mag
self.Ipmu_phase = Ipmu_phase
class Zdata_init:
def __init__(self, meas):
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)
self.mstddev = numpy.zeros(nmeas)
""" Insert here per unit values of the grid for power and voltage """
S = 100*(10**6)
V = (11*(10**3))/math.sqrt(3)
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)
""" Write here the indexes of the nodes/branches where the measurements are"""
V_idx = numpy.array([1,11,55])
I_idx = numpy.array([])
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([])
""" Write here the percent uncertainties of the measurements"""
V_unc = 1
I_unc = 0
Sinj_unc = 5
S_unc = 2
Pmu_mag_unc = 0.7
Pmu_phase_unc = 0.7
V = Measurements(V_idx,V_unc)
I = Measurements(I_idx,I_unc)
Sinj = Measurements(Sinj_idx,Sinj_unc)
S1 = Measurements(S1_idx,S_unc)
S2 = Measurements(S2_idx,S_unc)
Ipmu_mag = Measurements(Ipmu_idx,Pmu_mag_unc)
Ipmu_phase = Measurements(Ipmu_idx,Pmu_phase_unc)
Vpmu_mag = Measurements(Vpmu_idx,Pmu_mag_unc)
Vpmu_phase = Measurements(Vpmu_idx,Pmu_phase_unc)
meas = Measurement_set(V, I, Sinj, S1, S2, Ipmu_mag, Ipmu_phase, Vpmu_mag, Vpmu_phase)
zdata = Zdata_init(meas)
zdatameas = Zdata_init(meas)
zdata, zdatameas = Meas_data.Zdatatrue_creation(zdata, zdatameas, meas, branch, Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue)
MC_trials = 1
#Iest_matrix = numpy.zeros((MC_trials,branch.num),dtype=numpy.complex_)
#Vest_matrix = numpy.zeros((MC_trials,node.num),dtype=numpy.complex_)
#S1est_matrix = numpy.zeros((MC_trials,branch.num),dtype=numpy.complex_)
#S2est_matrix = numpy.zeros((MC_trials,branch.num),dtype=numpy.complex_)
#Sinjest_matrix = numpy.zeros((MC_trials,node.num),dtype=numpy.complex_)
iter_counter = numpy.zeros(MC_trials)
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)
# Iest_matrix[mciter] = Iest.complex
# Vest_matrix[mciter] = Vest.complex
# S1est_matrix[mciter] = S1est.complex
# S2est_matrix[mciter] = S2est.complex
# Sinjest_matrix[mciter] = Sinjest.complex
\ No newline at end of file
def Zdatatrue_creation(zdata, zdatameas, meas, branch, V, I, Iinj, S1, S2, Sinj):
""" It gives the true values at the measurement points."""
import numpy
t1 = numpy.ones(meas.V.num)
t2 = 2*numpy.ones(meas.Sinj.num)
t3 = 3*numpy.ones(meas.Sinj.num)
t4 = 4*numpy.ones(meas.S1.num + meas.S2.num)
t5 = 5*numpy.ones(meas.S1.num + meas.S2.num)
t6 = 6*numpy.ones(meas.I.num)
t7 = 7*numpy.ones(meas.Vpmu_mag.num)
t8 = 8*numpy.ones(meas.Vpmu_phase.num)
t9 = 9*numpy.ones(meas.Ipmu_mag.num)
t10 = 10*numpy.ones(meas.Ipmu_phase.num)
zdata.mtype = numpy.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
zdatameas.mtype = numpy.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
z1 = V.mag[meas.V.index-1]
z2 = Sinj.real[meas.Sinj.index-1]
z3 = Sinj.imag[meas.Sinj.index-1]
z4_1 = S1.real[meas.S1.index-1]
z4_2 = S2.real[meas.S2.index-1]
z5_1 = S1.imag[meas.S1.index-1]
z5_2 = S2.imag[meas.S2.index-1]
z6 = I.mag[meas.I.index-1]
z7 = V.mag[meas.Vpmu_mag.index-1]
z8 = V.phase[meas.Vpmu_phase.index-1]
z9 = I.mag[meas.Ipmu_mag.index-1]
z10 = I.phase[meas.Ipmu_phase.index-1]
zdata.mval = numpy.concatenate((z1,z2,z3,z4_1,z4_2,z5_1,z5_2,z6,z7,z8,z9,z10),axis=0)
b1 = numpy.zeros(meas.V.num)
b2 = numpy.zeros(meas.Sinj.num)
b3 = numpy.zeros(meas.Sinj.num)
b4_1 = branch.code[meas.S1.index-1]
b4_2 = -branch.code[meas.S2.index-1]
b5_1 = branch.code[meas.S1.index-1]
b5_2 = -branch.code[meas.S2.index-1]
b6 = branch.code[meas.I.index-1]
b7 = numpy.zeros(meas.Vpmu_mag.num)
b8 = numpy.zeros(meas.Vpmu_phase.num)
b9 = branch.code[meas.Ipmu_mag.index-1]
b10 = branch.code[meas.Ipmu_phase.index-1]
zdata.mbranch = numpy.concatenate((b1,b2,b3,b4_1,b4_2,b5_1,b5_2,b6,b7,b8,b9,b10),axis=0)
zdata.mbranch = zdata.mbranch.astype(int)
zdatameas.mbranch = numpy.concatenate((b1,b2,b3,b4_1,b4_2,b5_1,b5_2,b6,b7,b8,b9,b10),axis=0)
zdatameas.mbranch = zdata.mbranch.astype(int)
fr1 = meas.V.index
fr2 = meas.Sinj.index
fr3 = meas.Sinj.index
fr4_1 = branch.start[meas.S1.index-1]
fr4_2 = branch.end[meas.S2.index-1]
fr5_1 = branch.start[meas.S1.index-1]
fr5_2 = branch.end[meas.S2.index-1]
fr6 = branch.start[meas.I.index-1]
fr7 = meas.Vpmu_mag.index
fr8 = meas.Vpmu_phase.index
fr9 = branch.start[meas.Ipmu_mag.index-1]
fr10 = branch.start[meas.Ipmu_phase.index-1]
zdata.mfrom = numpy.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
zdata.mfrom = zdata.mfrom.astype(int)
zdatameas.mfrom = numpy.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
zdatameas.mfrom = zdata.mfrom.astype(int)
to1 = numpy.zeros(meas.V.num)
to2 = numpy.zeros(meas.Sinj.num)
to3 = numpy.zeros(meas.Sinj.num)
to4_1 = branch.end[meas.S1.index-1]
to4_2 = branch.start[meas.S2.index-1]
to5_1 = branch.end[meas.S1.index-1]
to5_2 = branch.start[meas.S2.index-1]
to6 = branch.end[meas.I.index-1]
to7 = numpy.zeros(meas.Vpmu_mag.num)
to8 = numpy.zeros(meas.Vpmu_phase.num)
to9 = branch.end[meas.Ipmu_mag.index-1]
to10 = branch.end[meas.Ipmu_phase.index-1]
zdata.mto = numpy.concatenate((to1,to2,to3,to4_1,to4_2,to5_1,to5_2,to6,to7,to8,to9,to10),axis=0)
zdata.mto = zdata.mto.astype(int)
zdatameas.mto = numpy.concatenate((to1,to2,to3,to4_1,to4_2,to5_1,to5_2,to6,to7,to8,to9,to10),axis=0)
zdatameas.mto = zdata.mto.astype(int)
d1 = z1*(meas.V.unc/300)
d2 = numpy.absolute(z2*(meas.Sinj.unc/300))
d3 = numpy.absolute(z3*(meas.Sinj.unc/300))
d4_1 = numpy.absolute(z4_1*(meas.S1.unc/300))
d4_2 = numpy.absolute(z4_2*(meas.S2.unc/300))
d5_1 = numpy.absolute(z5_1*(meas.S1.unc/300))
d5_2 = numpy.absolute(z5_2*(meas.S2.unc/300))
d6 = z6*(meas.I.unc/300)
d7 = z7*(meas.Vpmu_mag.unc/300)
d8 = (meas.Vpmu_phase.unc/300)*numpy.ones(meas.Vpmu_phase.num)
d9 = z9*(meas.Ipmu_mag.unc/300)
d10 = (meas.Ipmu_phase.unc/300)*numpy.ones(meas.Ipmu_phase.num)
zdata.mstddev = numpy.concatenate((d1,d2,d3,d4_1,d4_2,d5_1,d5_2,d6,d7,d8,d9,d10),axis=0)
zdatameas.mstddev = numpy.concatenate((d1,d2,d3,d4_1,d4_2,d5_1,d5_2,d6,d7,d8,d9,d10),axis=0)
# idx = numpy.where(zdata.mstddev<10**(-6))
# zdata.mstddev[idx] = 10**(-6)
return zdata, zdatameas
def Zdatameas_creation(zdata, zdatameas):
""" It gives the measured values (affected by uncertainty) at the measurement points."""
import numpy
zmeas = zdata.mval
zdev = zdata.mstddev
err_pu = numpy.random.normal(0,1,len(zmeas))
zdatameas.mval = zmeas + numpy.multiply(zdev,err_pu)
return zdatameas
\ No newline at end of file
def Network_95_nodes(Base, slackV):
import numpy
class Branch:
def __init__(self, bstart, bend, blength, br, bx, Base):
self.num = len(bstart)
self.start = bstart
self.end = bend
self.code = numpy.linspace(1,self.num,self.num)
self.R = numpy.multiply(br,blength)
self.X = numpy.multiply(bx,blength)
# self.R = numpy.divide(R,Base.Z)
# self.X = numpy.divide(X,Base.Z)
self.Z = self.R + 1j*self.X
class Node:
def __init__(self, ntype, slackV, P, Q):
self.type = ntype
self.P = P
self.Q = Q
P2 = numpy.concatenate(([slackV,],P),axis=0)
Q2 = numpy.concatenate(([0,],Q),axis=0)
self.pwr_flow_values = [ntype, P2, Q2]
self.num = len(ntype)
bstart = numpy.array([1, 2, 1, 4, 5, 5, 7, 7, 9, 9, 11, 12, 13, 13, 15, 15, 17, 17, 19, 20,
19, 22, 23, 24, 25, 25, 27, 11, 29, 29, 31, 31, 33, 34, 33, 36, 37, 38, 39, 37,
41, 42, 43, 43, 45, 46, 47, 48, 48, 50, 50, 52, 52, 42, 55, 56, 56, 58, 55, 60,
61, 62, 61, 64, 65, 65, 67, 68, 69, 69, 71, 72, 73, 74, 74, 60, 77, 78, 79, 80,
78, 82, 82, 84, 85, 86, 84, 88, 89, 90, 91, 92, 93, 94])
bend = numpy.array([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])
bR = numpy.array([0.05489, 0.03881, 0.04879, 0.09755, 0.17322, 0.21, 0.24251, 0.2586, 0.34645, 0.1293, 0.30169, 0.19395,
0.17322, 0.23705, 0.20787, 0.2586, 0.13858, 0.1293, 0.1724, 0.10775, 0.2155, 0.19395, 0.30169, 0.09,
0.13161, 0.36, 0.12, 0.295, 0.20787, 0.354, 0.24251, 0.27716, 0.2586, 0.3118, 0.11149, 0.34612, 0.15608,
0.22298, 0.2135, 0.34479, 0.1293, 0.118, 0.20787, 0.236, 0.177, 0.09401, 0.177, 0.236, 0.354, 0.354,
0.27716, 0.2135, 0.53374, 0.1724, 0.20787, 0.27716, 0.41574, 0.27716, 0.30169, 0.36517, 0.17322, 0.3118,
0.14607, 0.25562, 0.20787, 0.18258, 0.29213, 0.25562, 0.09401, 0.4382, 0.2191, 0.07521, 0.14607, 0.29213,
0.40168, 0.2586, 0.1293, 0.18258, 0.29213, 0.4382, 0.2155, 0.1293, 0.19395, 0.27716, 0.48502, 0.22621,
0.2586, 0.06465, 0.27244, 0.16346, 0.0862, 0.15085, 0.27244, 0.49039])
bX = numpy.array([0.0569, 0.104, 0.05058, 0.33284, 0.07589, 0.203, 0.10624, 0.17673, 0.15178, 0.08836, 0.20618, 0.13255,
0.07589, 0.172, 0.09107, 0.17673, 0.06071, 0.08836, 0.11782, 0.07364, 0.14727, 0.13255, 0.20618, 0.087,
0.05033, 0.348, 0.116, 0.15, 0.09107, 0.18, 0.10624, 0.12142, 0.17673, 0.1366, 0.07376, 0.20653, 0.10326,
0.14752, 0.09126, 0.23564, 0.08836, 0.06, 0.09107, 0.12, 0.09, 0.03595, 0.09, 0.12, 0.18, 0.18, 0.12142,
0.09126, 0.22816, 0.11782, 0.09107, 0.12142, 0.18213, 0.12142, 0.20618, 0.15244, 0.07589, 0.1366, 0.06098,
0.10671, 0.09107, 0.07622, 0.12195, 0.10671, 0.03595, 0.18293, 0.09146, 0.02876, 0.06098, 0.12195, 0.16768,
0.17673, 0.08836, 0.07622, 0.12195, 0.17673, 0.14727, 0.08836, 0.13255, 0.12142, 0.21249, 0.04686, 0.17673,
0.04418, 0.04012, 0.02407, 0.05891, 0.10309, 0.04012, 0.07222])
blength = numpy.linspace(1,1,94)
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,
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.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))
ntype = ['slack']
for j in range(1, len(P)+1):
ntype.append(nPQ)
branch = Branch(bstart, bend, blength, bR, bX, Base)
node = Node(ntype, slackV, P, Q)
return branch, node
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
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 BcDsseCall(branch, node, zdata):
""" It identifies the type of measurements present in the measurement set and
calls the appropriate estimator for dealing with them."""
trad_code = 0
PMU_code = 0
Vmag_meas = numpy.where(zdata.mtype==1)
if len(Vmag_meas[0])>0:
trad_code = 1
Vpmu_meas = numpy.where(zdata.mtype==7)
if len(Vpmu_meas[0])>0:
PMU_code = 2
est_code = trad_code + PMU_code
if est_code == 1:
Vest, Iest = BcDsseTrad(branch, node, zdata)
elif est_code == 2:
Vest, Iest = BcDssePmu(branch, node, zdata)
else:
Vest, Iest = BcDsseMixed(branch, node, zdata)
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(Iest.real[to[0]]) - numpy.sum(Iest.real[fr[0]])
Iinj_x[k-1] = numpy.sum(Iest.imag[to[0]]) - numpy.sum(Iest.imag[fr[0]])
Iinjest = Real_to_all(Iinj_r,Iinj_x)
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))
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):
""" It performs state estimation using branch current state variables and it is customized
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)
pfidx = numpy.where(zdata.mtype==4)
qfidx = numpy.where(zdata.mtype==5)
iidx = numpy.where(zdata.mtype==6)
nvi = len(vidx[0])
npi = len(pidx[0])
nqi = len(qidx[0])
npf = len(pfidx[0])
nqf = len(qfidx[0])
nii = len(iidx[0])
busvi = zdata.mfrom[vidx]
buspi = zdata.mfrom[pidx]
busqi = zdata.mfrom[qidx]
buspf = zdata.mfrom[pfidx]
busqf = zdata.mfrom[qfidx]
branchpf = zdata.mbranch[pfidx]
branchqf = zdata.mbranch[qfidx]
branchii = zdata.mbranch[iidx]
z = zdata.mval
idx = numpy.where(zdata.mbranch<0)
z[idx] = -z[idx]
Pinj = z[pidx]