Commit 6161a224 authored by Jan Dinkelbach's avatar Jan Dinkelbach
Browse files

moved python version of marco pau's original implementation to separate branch

parent 0a0265c5
# State Estimation Algorithm
The algorithm is implemented in Python and derived from Marco Pau's Matlab wonderful version.
The algorithm is implemented in Python and derived from Marco Pau's original version.
This project uses CIMpy to read CIM data into Python objects: https://git.rwth-aachen.de/acs/core/cim/cimpy
\ No newline at end of file
__all__ = ["bc_powerflow", "bc_state_estimator", "nv_powerflow", "nv_state_estimator", "network", "nv_powerflow_cim", "results"]
\ 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)
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)
self.R = self.Z.real
self.X = self.Z.imag
def solve(branches, nodes):
"""It performs Power Flow by using rectangular branch current state variables."""
znod = Znod(branches,nodes)
# real matrices storing complex numbers
# real in position n and imag in position n+1
z = np.zeros(2*(branches.num+1))
h = np.zeros(2*(branches.num+1))
H = np.zeros((2*(branches.num+1),2*(branches.num+1)))
for k in range(1,nodes.num+1):
i = k-1
m = 2*i
if nodes.type[i] == 'slack':
V = nodes.pwr_flow_values[1][i]
theta = nodes.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:] = 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 nodes.type[i] == 'PQ':
to = np.where(branches.end==k)
fr = np.where(branches.start==k)
to1 = to[0]+2
to2 = to1+branches.num
fr1 = fr[0]+2
fr2 = fr1+branches.num
H[m][to1] = 1
H[m+1][to2] = 1
H[m][fr1] = -1
H[m+1][fr2] = -1
elif nodes.type[i] == 'PV':
z[m+1] = nodes.pwr_flow_values[2][i]
to = np.where(branches.end==k)
fr = np.where(branches.start==k)
to1 = to[0]+2
fr1 = fr[0]+2
H[m][to1] = 1
H[m][fr1] = -1
diff = 5
epsilon = 10**(-10)
V = np.ones(nodes.num) + 1j* np.zeros(nodes.num)
num_iter = 0
State = np.ones(2*branches.num)
State = np.concatenate((np.array([1,0]),State),axis=0)
while diff > epsilon:
for k in range(1,nodes.num+1):
i = k-1
m = 2*i
if nodes.type[i] == 'slack':
h[m] = np.inner(H[m],State)
h[m+1] = np.inner(H[m+1],State)
elif nodes.type[i] == 'PQ':
z[m] = (nodes.pwr_flow_values[1][i]*V[i].real + nodes.pwr_flow_values[2][i]*V[i].imag)/(np.abs(V[i])**2)
z[m+1] = (nodes.pwr_flow_values[1][i]*V[i].imag - nodes.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 nodes.type[i] == 'PV':
z[m] = (nodes.pwr_flow_values[1][i]*V[i].real + nodes.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]+branches.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
diff = np.amax(np.absolute(Delta_State))
V[0] = State[0] + 1j * State[1]
Ir = State[2:branches.num+2]
Ix = State[branches.num+2:]
Irx = Ir + 1j*Ix
Vcomplex = np.inner((V[0].real + 1j*V[0].imag), np.ones(nodes.num))
DeltaV = np.inner(znod.Z,Irx)
V = Vcomplex[0] - DeltaV
num_iter = num_iter+1
I = Ir + 1j*Ix
Iinj_r = np.zeros(nodes.num)
Iinj_x = np.zeros(nodes.num)
for k in range(1,nodes.num+1):
to = np.where(branches.end==k)
fr = np.where(branches.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 = np.multiply(V[branches.start-1], np.conj(I))
S2 = np.multiply(V[branches.end-1], np.conj(I))
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):
""" 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 DsseCall(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 = DsseTrad(branch, node, zdata)
elif est_code == 2:
Vest, Iest = DssePmu(branch, node, zdata)
else:
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):
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 DsseTrad(branch, node, zdata):
""" It performs state estimation using branch current state variables and it is customized
to work without PMU measurements. """
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]
Qinj = z[qidx]
Pbr = z[pfidx]
Qbr = z[qfidx]
idx = numpy.where(zdata.mstddev<10**(-6))
zdata.mstddev[idx] = 10**(-6)
weights = zdata.mstddev**(-2)
W = numpy.diag(weights)
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):
m = buspi[i]
to = numpy.where(branch.end==m)
fr = numpy.where(branch.start==m)
to1 = to[0]+1
to2 = to1+branch.num
fr1 = fr[0]+1
fr2 = fr1+branch.num
H2[i][to1] = 1
H2[i][fr1] = -1
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):
br = abs(branchpf[i])
H4[i][br] = 1
H5[i][br+branch.num] = 1
epsilon = 5
Vr = numpy.ones(node.num)
Vx = numpy.zeros(node.num)
V = Real_to_all(Vr, Vx)
Ir = numpy.zeros(node.num)
Ix = numpy.zeros(node.num)
I = Real_to_all(Ir, Ix)
num_iter = 0
State = numpy.zeros(2*branch.num)
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
z[qidx] = Ixinj
Irbr = (Pbr*V.real[buspf-1] + Qbr*V.imag[busqf-1])/(V.mag[buspf-1]**2)
Ixbr = (Pbr*V.imag[buspf-1] - Qbr*V.real[busqf-1])/(V.mag[buspf-1]**2)
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):
m = busvi[i]-1
if m == 0:
H1[i][0] = 1
else:
idx = numpy.where(znod.Z[m]!=0)
H1[i][0] = numpy.cos(V.phase[m])
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))
for i in range(0,nii):
idx = branchii[i]
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)
g = numpy.inner(H.transpose(),numpy.inner(W,res))
WH = numpy.inner(W,H.transpose())
G = numpy.inner(H.transpose(),WH.transpose())
Ginv = numpy.linalg.inv(G)
Delta_State = numpy.inner(Ginv,g)
State = State + Delta_State
epsilon = numpy.amax(numpy.absolute(Delta_State))
V.real[0] = State[0]
Ir = State[1:branch.num+1]
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
V.real = numpy.real(V.complex)
V.imag = numpy.imag(V.complex)
V = Real_to_all(V.real, V.imag)
I = Real_to_all(Ir,Ix)
num_iter = num_iter+1
return V, I
def DssePmu(branch, node, zdata):
""" It performs state estimation using branch current state variables and it is customized
to work using only 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
znod = Znod(branch,node)
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."""
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
znod = Znod(branch,node)
\ No newline at end of file
import numpy as np
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 = np.zeros(nmeas)
self.mval = np.zeros(nmeas)
self.mbranch = np.zeros(nmeas)
self.mfrom = np.zeros(nmeas)
self.mto = np.zeros(nmeas)
self.mstddev = np.zeros(nmeas)
def Zdatatrue_creation(zdata, zdatameas, meas, system, V, I, Iinj, S1, S2, Sinj):
""" It gives the true values at the measurement points."""
import numpy as np
t1 = np.ones(meas.V.num)
t2 = 2*np.ones(meas.Sinj.num)
t3 = 3*np.ones(meas.Sinj.num)
t4 = 4*np.ones(meas.S1.num + meas.S2.num)
t5 = 5*np.ones(meas.S1.num + meas.S2.num)
t6 = 6*np.ones(meas.I.num)
t7 = 7*np.ones(meas.Vpmu_mag.num)
t8 = 8*np.ones(meas.Vpmu_phase.num)
t9 = 9*np.ones(meas.Ipmu_mag.num)
t10 = 10*np.ones(meas.Ipmu_phase.num)
zdata.mtype = np.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
zdatameas.mtype = np.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
z1 = np.abs(V[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 = np.abs(I[meas.I.index-1])
z7 = np.abs(V[meas.Vpmu_mag.index-1])
z8 = np.angle(V[meas.Vpmu_phase.index-1])
z9 = np.abs(I[meas.Ipmu_mag.index-1])
z10 = np.angle(I[meas.Ipmu_phase.index-1])
zdata.mval = np.concatenate((z1,z2,z3,z4_1,z4_2,z5_1,z5_2,z6,z7,z8,z9,z10),axis=0)
"""
b1 = np.zeros(meas.V.num)
b2 = np.zeros(meas.Sinj.num)
b3 = np.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 = np.zeros(meas.Vpmu_mag.num)
b8 = np.zeros(meas.Vpmu_phase.num)
b9 = branch.code[meas.Ipmu_mag.index-1]
b10 = branch.code[meas.Ipmu_phase.index-1]
zdata.mbranch = np.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 = np.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 = np.array([])
fr4_2 = np.array([])
fr5_1 = np.array([])
fr5_2 = np.array([])
for index in range(len(meas.S1.index)):
fr4_1 = np.append(fr4_1, system.branches[meas.S1.index[index]-1].start_node.index)
fr5_1 = np.append(fr5_1, system.branches[meas.S1.index[index]-1].start_node.index)
for index in range(len(meas.S2.index)):
fr4_2 = np.append(fr4_2, system.branches[meas.S2.index[index]-1].end_node.index)
fr5_2 = np.append(fr5_2, system.branches[meas.S2.index[index]-1].end_node.index)
fr6 = np.array([])
for index in range(len(meas.I.index)):
fr6 = np.append(fr6, system.branches[meas.I.index[index]-1].start_node.index)
fr7 = meas.Vpmu_mag.index
fr8 = meas.Vpmu_phase.index
fr9 = np.array([])
fr10 = np.array([])
for index in range(len(meas.Ipmu_mag.index)):
fr9 = np.append(fr9, system.branches[meas.Ipmu_mag.index[index]-1].start_node.index)
for index in range(len(meas.Ipmu_phase.index)):
fr10 = np.append(fr10, system.branches[meas.Ipmu_phase.index[index]-1].start_node.index)
zdata.mfrom = np.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 = np.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 = np.zeros(meas.V.num)
to2 = np.zeros(meas.Sinj.num)
to3 = np.zeros(meas.Sinj.num)
to4_1 = np.array([])
to4_2 = np.array([])
to5_1 = np.array([])
to5_2 = np.array([])
for index in range(len(meas.S1.index)):
to4_1 = np.append(to4_1, system.branches[meas.S1.index[index]-1].end_node.index)
to5_1 = np.append(to5_1, system.branches[meas.S1.index[index]-1].end_node.index)