Skip to content
GitLab
Commit 5cb610b2 authored by Markus Mirz's avatar Markus Mirz
Browse files

add node voltage powerflow and estimator

parent f62b5eb3
Expand all Show whitespace changes
Inline Side-by-side
acs-state-estimation/Main.py
View file @ 5cb610b2
......@@ -2,9 +2,10 @@ import numpy
import math
import matplotlib as plt
import Network_data
import powerflow
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 = powerflow.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
......
acs-state-estimation/Network_data.py
View file @ 5cb610b2
......@@ -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):
......
acs-state-estimation/StateEstimator.py → acs-state-estimation/bc_state_estimator.py
View file @ 5cb610b2
......@@ -2,13 +2,38 @@ 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)
def BcDsseCall(branch, node, zdata):
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."""
......@@ -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,33 +70,16 @@ 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"""
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
to work without PMU measurements. """
vidx = numpy.where(zdata.mtype==1)
pidx = numpy.where(zdata.mtype==2)
......@@ -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):
......
acs-state-estimation/nv_powerflow.py 0 → 100644
View file @ 5cb610b2
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 rectangular branch current state variables."""
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)
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.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.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)
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
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
acs-state-estimation/nv_state_estimator.py 0 → 100644
View file @ 5cb610b2
This diff is collapsed.