Commit 502c02b9 authored by martin.moraga's avatar martin.moraga
Browse files

improve nv_state_estimator_cim.py

parent 465864e2
......@@ -7,26 +7,6 @@ from network import BusType
sys.path.append("../../../dataprocessing")
from villas.dataprocessing.readtools import read_timeseries_dpsim
class PF_Results():
def __init__(self):
self.V = [] #Node Voltage
self.I = [] #Branch Current
self.Iinj = [] #Injection Current
self.S1 = [] #Branch Complex Power measured at 1st node of the line
self.S2 = [] #Branch Complex Power measured at 2nd node of the line
self.SInj = [] #Complex Power Injection
self.num_iter = 0 #Number of Iterations of Newton Rapson
def read_data(self, file_name, system):
"""
To read the voltages from a input file
"""
loadflow_results = read_data(loadflow_results_file)
#order readed data according to system.nodes
self.V = np.zeros(len(loadflow_results), dtype=np.complex_)
for elem in range(len(system.nodes)):
self.V[elem] = loadflow_results[system.nodes[elem].uuid]
class PowerflowNode():
def __init__(self, topo_node):
self.topology_node = topo_node
......@@ -38,8 +18,9 @@ class PowerflowBranch():
def __init__(self, topo_branch):
self.topology_branch = topo_branch
self.current = complex(0, 0)
self.power = complex(0, 0)
self.power = complex(0, 0) #complex power flow at branch, measured at initial node
self.power2 = complex(0, 0) #complex power flow at branch, measured at final node
class PowerflowResults():
def __init__(self, system):
self.nodes=[]
......@@ -51,13 +32,13 @@ class PowerflowResults():
for branch in system.branches:
self.branches.append(PowerflowBranch(topo_branch=branch))
def read_data(self, file_name):
def read_data_dpsim(self, file_name):
"""
read the voltages from a input file
read the voltages from a dpsim input file
"""
loadflow_results = read_data(loadflow_results_file)
for elem in range(len(self.nodes)):
elem.V = loadflow_results[elem.topology_node.uuid]
loadflow_results = read_timeseries_dpsim(file_name, print_status=False)
for node in self.nodes:
node.V = loadflow_results[node.topology_node.uuid].values[0]
def load_voltages(self, V):
"""
......@@ -98,9 +79,29 @@ class PowerflowResults():
for node in self.nodes:
node.power = node.voltage*np.conj(node.current)
def calculateS1(self):
"""
calculate complex power flow at branch, measured at initial node
"""
for branch in self.branches:
branch_index = branch.topology_branch.start_node.index
for node in self.nodes:
if branch_index == node.topology_node.index:
branch.power = node.voltage*(np.conj(branch.current))
def calculateS2(self):
"""
calculate complex ower flow at branch, measured at final node
"""
for branch in self.branches:
branch_index = branch.topology_branch.end_node.index
for node in self.nodes:
if branch_index == node.topology_node.index:
branch.power2 = -node.voltage*(np.conj(branch.current))
def get_voltages(self):
"""
get node voltages
get complex Power Injection at nodes
for a test purpose
"""
voltages = np.zeros(len(self.nodes), dtype=np.complex_)
......@@ -134,16 +135,34 @@ class PowerflowResults():
get branch currents
for a test purpose
"""
currents = np.zeros(len(self.branches), dtype=np.complex_)
for elem in self.branches:
voltages[elem.topology_node.index] = elem.voltage
return voltages
I = np.array([])
for branch in self.branches:
I = np.append(I, branch.current)
return I
def get_S1(self):
"""
get complex Power flow at branch, measured at initial node
for a test purpose
"""
S1 = np.array([])
for branch in self.branches:
S1 = np.append(S, branch.power)
return S1
def get_S2(self):
"""
get complex Power flow at branch, measured at final node
for a test purpose
"""
S2 = np.array([])
for branch in self.branches:
S2 = np.append(S, branch.power2)
return S2
def solve(system):
"""It performs Power Flow by using rectangular node voltage state variables."""
Ymatrix, Adj = Ymatrix_calc(system)
nodes_num = len(system.nodes)
branches_num = len(system.branches)
......@@ -161,24 +180,24 @@ def solve(system):
H[m][i] = 1
H[m+1][i2] = 1
elif type is BusType.PQ:
H[m][i] = - np.real(Ymatrix[i][i])
H[m][i2] = np.imag(Ymatrix[i][i])
H[m+1][i] = - np.imag(Ymatrix[i][i])
H[m+1][i2] = - np.real(Ymatrix[i][i])
idx1 = np.subtract(Adj[i],1)
H[m][i] = - np.real(system.Ymatrix[i][i])
H[m][i2] = np.imag(system.Ymatrix[i][i])
H[m+1][i] = - np.imag(system.Ymatrix[i][i])
H[m+1][i2] = - np.real(system.Ymatrix[i][i])
idx1 = np.subtract(system.Adjacencies[i],1)
idx2 = idx1 + nodes_num
H[m][idx1] = - np.real(Ymatrix[i][idx1])
H[m][idx2] = np.imag(Ymatrix[i][idx1])
H[m+1][idx1] = - np.imag(Ymatrix[i][idx1])
H[m+1][idx2] = - np.real(Ymatrix[i][idx1])
H[m][idx1] = - np.real(system.Ymatrix[i][idx1])
H[m][idx2] = np.imag(system.Ymatrix[i][idx1])
H[m+1][idx1] = - np.imag(system.Ymatrix[i][idx1])
H[m+1][idx2] = - np.real(system.Ymatrix[i][idx1])
elif type is BusType.PV:
z[m+1] = np.real(system.nodes[i].power)
H[m][i] = - np.real(Ymatrix[i][i])
H[m][i2] = np.imag(Ymatrix[i][i])
idx1 = np.subtract(Adj[i],1)
H[m][i] = - np.real(system.Ymatrix[i][i])
H[m][i2] = np.imag(system.Ymatrix[i][i])
idx1 = np.subtract(system.Adjacencies[i],1)
idx2 = idx1 + nodes_num
H[m][idx1] = - np.real(Ymatrix[i][idx1])
H[m][idx2] = np.imag(Ymatrix[i][idx1])
H[m][idx1] = - np.real(system.Ymatrix[i][idx1])
H[m][idx2] = np.imag(system.Ymatrix[i][idx1])
epsilon = 10**(-10)
diff = 5
......@@ -215,10 +234,18 @@ def solve(system):
diff = np.amax(np.absolute(Delta_State))
V = State[:nodes_num] + 1j * State[nodes_num:]
num_iter = num_iter+1
return V, num_iter
results = PowerflowResults(system)
results.load_voltages(V)
results.calculateI()
results.calculateIinj()
results.calculateSinj()
results.calculateI()
results.calculateS()
return results, num_iter
def calculateI(system, V):
"""
......
import sys
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 Complex_to_all:
def __init__(self,Arx):
""" Converts complex numbers to all the other formats. """
self.complex = Arx
self.real = numpy.real(self.complex)
self.imag = numpy.imag(self.complex)
self.mag = numpy.absolute(self.complex)
self.phase = numpy.angle(self.complex)
sys.path.append("../../../dataprocessing")
from villas.dataprocessing.readtools import read_timeseries_dpsim
class Measurement():
def __init__(self, topo_node, meas_type, meas_value, std_dev):
"""
Creates a measurement, which is used by the estimation module. Possible types of measurements are: v, p, q, i, Vpmu and Ipmu
@meas_type: "p" for active power, "q" for reactive power, "v" for voltage, "i" for electrical current, "Vpmu" and "Ipmu"
@node_uuid:
@meas_value: measurement value.
@std_dev: standard deviation in the same unit as the measurement.
"""
if meas_type not in ("v", "Vpmu", "p", "q", "i", "Ipmu"):
raise Exception("Invalid measurement type (%s)" % meas_type)
def DsseCall(system, zdata, Ymatrix, Adj):
self.topology_node = topo_node
self.meas_type = meas_type
self.meas_value = meas_value
self.std_dev = self.meas_value*(std_dev/300)
class Measurents_set():
def __init__(self):
self.measurements_set = []
def create_measurement(self, topo_node, meas_type, meas_value, std_dev):
self.measurements_set.append(Measurement(topo_node, meas_type, meas_value, std_dev))
def DsseCall(system, zdata, measurements):
""" It identifies the type of measurements present in the measurement set and
calls the appropriate estimator for dealing with them."""
# system: model of the system (nodes, lines, topology)
# zdata: Vector of measurements in Input (voltages, currents, powers)
# Ymatrix: Parameters of the lines (R,X,G,B)
# Adj: Adjacency Matrix
# 1. Select type of Estimator.
# If at least a PMU is present we launch the combined estimator, otherwise a simple traditional etimator
trad_code = 0
PMU_code = 0
Vmag_meas = numpy.where(zdata.mtype==1) #Voltage Magnitude measurement: type 1
if len(Vmag_meas[0])>0:
trad_code = 1
Vpmu_meas = numpy.where(zdata.mtype==7) #Voltage Phasor (Magnitude+Phase) measurement: type 7
if len(Vpmu_meas[0])>0:
PMU_code = 2
Vmag_meas = 0
Vpmu_meas = 0
for elem in measurements.measurements_set:
if elem.meas_type=="v":
Vmag_meas = Vmag_meas + 1
elif elem.meas_type=="Vpmu":
Vpmu_meas = Vpmu_meas + 1
trad_code = 1 if Vmag_meas >0 else 0
PMU_code = 2 if Vpmu_meas >0 else 0
print(trad_code)
print(PMU_code)
est_code = trad_code + PMU_code
# 2. Run Estimator.
if est_code == 1:
Vest = DsseTrad(system, zdata, Ymatrix, Adj)
Vest = DsseTrad(system, zdata)
elif est_code == 2:
Vest = DssePmu(system, zdata, Ymatrix, Adj)
Vest = DssePmu(system, zdata)
else:
Vest = DsseMixed(system, zdata, Ymatrix, Adj)
Vest = DsseMixed(system, zdata)
# 3. Populate arrays with data in output of estimator
# The Estimator provides:
......@@ -64,15 +74,15 @@ def DsseCall(system, zdata, Ymatrix, Adj):
branches_num = len(system.branches)
""" From here on, all the other quantities of the grid are calculated """
Irx = numpy.zeros(branches_num, dtype=numpy.complex)
Iest = numpy.zeros(branches_num, dtype=numpy.complex)
for index in range(branches_num):
fr = system.branches[index].start_node.index #branch.start[idx]-1
to = system.branches[index].end_node.index #branch.end[idx]-1
Irx[index] = - (Vest.complex[fr] - Vest.complex[to])*Ymatrix[fr][to]
Irx[index] = - (Vest.complex[fr] - Vest.complex[to])*system.Ymatrix[fr][to]
Ir = numpy.real(Irx)
Ix = numpy.imag(Irx)
Iest = Real_to_all(Ir,Ix)
#Iest = Real_to_all(Ir,Ix)
Iinj_r = numpy.zeros(nodes_num)
Iinj_x = numpy.zeros(nodes_num)
for k in range(nodes_num):
......@@ -84,35 +94,31 @@ def DsseCall(system, zdata, Ymatrix, Adj):
if k==system.branches[m].end_node.index:
to.append(m)
Iinj_r[k] = numpy.sum(Iest.real[to]) - numpy.sum(Iest.real[fr])
Iinj_x[k] = numpy.sum(Iest.imag[to]) - numpy.sum(Iest.imag[fr])
Iinj_r[k] = numpy.sum(Iest[to].real) - numpy.sum(Iest[fr].real)
Iinj_x[k] = numpy.sum(Iest[to].imag) - numpy.sum(Iest[fr].imag)
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))
#Iinjest = Real_to_all(Iinj_r,Iinj_x)
Iinjest = Iinj_r + 1j*Iinj_x
Sinjest = numpy.multiply(Vest.complex,numpy.conj(Iinjest.complex))
#Sinjest = Real_to_all(numpy.real(Sinj_rx),numpy.imag(Sinj_rx))
S1_rx=numpy.array([])
S2_rx=numpy.array([])
S1est=numpy.array([])
S2est=numpy.array([])
for i in range(branches_num):
S1_rx=numpy.append(S1_rx, Vest.complex[system.branches[i].start_node.index]*numpy.conj(Iest.complex[i]))
S2_rx=numpy.append(S2_rx, -Vest.complex[system.branches[i].end_node.index]*numpy.conj(Iest.complex[i]))
S1est=numpy.append(S1_rx, Vest.complex[system.branches[i].start_node.index]*numpy.conj(Iest[i]))
S2est=numpy.append(S2_rx, -Vest.complex[system.branches[i].end_node.index]*numpy.conj(Iest[i]))
S1est = Real_to_all(numpy.real(S1_rx),numpy.imag(S1_rx))
S2est = Real_to_all(numpy.real(S2_rx),numpy.imag(S2_rx))
#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(system, zdata, Ymatrix, Adj):
def DsseTrad(system, zdata):
""" It performs state estimation using rectangular node voltage state variables
and it is customized to work without PMU measurements"""
# Traditional state estimator
# system: model of the system (nodes, lines, topology)
# zdata: Vector of measurements in Input (voltages, currents, powers)
# Ymatrix: Parameters of the lines (R,X,G,B)
# Adj: Adjacency Matrix
# zdata: Vector of measurements in Input (voltages, currents, powers)
nodes_num = len(system.nodes) # number of nodes of the grids, identify also the number of states (2*nodes_num-1)
......@@ -158,11 +164,10 @@ def DsseTrad(system, zdata, Ymatrix, Adj):
W = numpy.diag(weights)
# Admittances of the lines of the network
Admittance = Complex_to_all(Ymatrix)
Gmatrix = Admittance.real
Bmatrix = Admittance.imag
Yabs_matrix = Admittance.mag
Yphase_matrix = Admittance.phase
Gmatrix = system.Ymatrix.real
Bmatrix = system.Ymatrix.imag
Yabs_matrix = numpy.absolute(system.Ymatrix)
Yphase_matrix = numpy.angle(system.Ymatrix)
# Jacobian Matrix. Includes the derivatives of the measurements (voltages, currents, powers) with respect to the states (voltages)
......@@ -178,7 +183,7 @@ def DsseTrad(system, zdata, Ymatrix, Adj):
H2[i][m2] = Bmatrix[m][m]
H3[i][m] = - Bmatrix[m][m]
H3[i][m2] = - Gmatrix[m][m]
idx = numpy.subtract(Adj[m],1)
idx = numpy.subtract(system.Adj[m],1)
H2[i][idx] = - Gmatrix[m][idx]
H3[i][idx] = - Bmatrix[m][idx]
if 0 in idx:
......@@ -212,7 +217,7 @@ def DsseTrad(system, zdata, Ymatrix, Adj):
epsilon = 5 # treshold to stop Netwon Rapson iterations
Vr = numpy.ones(nodes_num) #initialize voltage real part to 1 per unit
Vx = numpy.zeros(nodes_num) #initialize voltage imaginary part to 0 per unit
V = Real_to_all(Vr, Vx)
V = Vr + 1j*Vx
num_iter = 0 # number of iterations of Newton Rapson
StateVr = numpy.ones(nodes_num) #initialize voltage real part to 1 per unit
......@@ -223,27 +228,27 @@ def DsseTrad(system, zdata, Ymatrix, Adj):
while epsilon>10**(-6):
""" Computation of equivalent current measurements in place of the power measurements """
# in every iteration the input power measurements are converted into currents by dividing by the voltage estimated at the previous iteration
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)
Irinj = (Pinj*V[buspi-1].real + Qinj*V[busqi-1].imag)/(numpy.absolute(V[buspi-1])**2)
Ixinj = (Pinj*V[buspi-1].imag - Qinj*V[busqi-1].real)/(numpy.absolute(V[buspi-1])**2)
z[pidx] = Irinj
z[qidx] = Ixinj
Irbr = (Pbr*V.real[fbuspf-1] + Qbr*V.imag[fbusqf-1])/(V.mag[fbuspf-1]**2)
Ixbr = (Pbr*V.imag[fbuspf-1] - Qbr*V.real[fbusqf-1])/(V.mag[fbuspf-1]**2)
Irbr = (Pbr*V[fbuspf-1].real + Qbr*V[fbusqf-1].imag)/(numpy.absolute(V[fbuspf-1])**2)
Ixbr = (Pbr*V[fbuspf-1].imag - Qbr*V[fbusqf-1].real)/(numpy.absolute(V[fbuspf-1])**2)
z[pfidx] = Irbr
z[qfidx] = Ixbr
""" Voltage Magnitude Measurements """
# at every iteration we update h(x) vector where V measure are available
h1 = V.mag[busvi-1]
h1 = V[busvi-1].mag[busvi-1]
# the Jacobian rows where voltage measurements are presents is updated
H1 = numpy.zeros((nvi,2*nodes_num-1))
for i in range(nvi):
m = busvi[i]-1
H1[i][m] = numpy.cos(V.phase[m])
H1[i][m] = numpy.cos(numpy.angle(V[m]))
if m > 0:
m2 = m + nodes_num-1
H1[i][m2] = numpy.sin(V.phase[m])
H1[i][m2] = numpy.sin(numpy.angle(V[m]))
""" Power Injection Measurements """
# h(x) vector where power injections are present
......@@ -265,8 +270,8 @@ def DsseTrad(system, zdata, Ymatrix, Adj):
for i in range(nii):
m = fbusiamp[i]-1
n = tbusiamp[i]-1
h6re[i] = Yabs_matrix[m][n]*((V.real[n]-V.real[m])*numpy.cos(Yphase_matrix[m][n]) + (V.imag[m]-V.imag[n])*numpy.sin(Yphase_matrix[m][n]))
h6im[i] = Yabs_matrix[m][n]*((V.real[n]-V.real[m])*numpy.sin(Yphase_matrix[m][n]) + (V.imag[n]-V.imag[m])*numpy.cos(Yphase_matrix[m][n]))
h6re[i] = Yabs_matrix[m][n]*((V[n].real-V[m].real)*numpy.cos(Yphase_matrix[m][n]) + (V[m].imag-V[n].imag)*numpy.sin(Yphase_matrix[m][n]))
h6im[i] = Yabs_matrix[m][n]*((V[n].real-V[m].real)*numpy.sin(Yphase_matrix[m][n]) + (V[n].imag-V[m].imag)*numpy.cos(Yphase_matrix[m][n]))
h6complex[i] = h6re[i] + 1j*h6im[i]
if num_iter>0:
h6[i] = numpy.absolute(h6complex[i])
......@@ -303,15 +308,13 @@ def DsseTrad(system, zdata, Ymatrix, Adj):
V.real = State[:nodes_num]
V.imag = State[nodes_num:]
V.imag = numpy.concatenate(([0],V.imag), axis=0)
V = Real_to_all(V.real, V.imag)
#V = Real_to_all(V.real, V.imag)
num_iter = num_iter+1
return V
def DssePmu(system, zdata, Ymatrix, Adj):
def DssePmu(system, zdata):
""" It performs state estimation using rectangular node voltage state variables
and it is customized to work using PMU measurements."""
......@@ -354,10 +357,9 @@ def DssePmu(system, zdata, Ymatrix, Adj):
weights = zdata.mstddev**(-2)
W = numpy.diag(weights)
Admittance = Complex_to_all(Ymatrix)
Gmatrix = Admittance.real
Bmatrix = Admittance.imag
Gmatrix = system.Ymatrix.real
Bmatrix = system.Ymatrix.imag
""" Jacobian for Power Injection Measurements (converted to equivalent
rectangualar current measurements) """
H2 = numpy.zeros((npi,2*nodes_num))
......@@ -369,7 +371,7 @@ def DssePmu(system, zdata, Ymatrix, Adj):
H2[i][m2] = Bmatrix[m][m]
H3[i][m] = - Bmatrix[m][m]
H3[i][m2] = - Gmatrix[m][m]
idx = numpy.subtract(Adj[m],1)
idx = numpy.subtract(system.Adj[m],1)
H2[i][idx] = - Gmatrix[m][idx]
H3[i][idx] = - Bmatrix[m][idx]
idx2 = idx + nodes_num
......@@ -449,7 +451,7 @@ def DssePmu(system, zdata, Ymatrix, Adj):
epsilon = 5
Vr = numpy.ones(nodes_num)
Vx = numpy.zeros(nodes_num)
V = Real_to_all(Vr, Vx)
V = Vr + 1j*Vx
num_iter = 0
StateVr = numpy.ones(nodes_num)
......@@ -461,16 +463,15 @@ def DssePmu(system, zdata, Ymatrix, Adj):
G = numpy.inner(H.transpose(),WH.transpose())
Ginv = numpy.linalg.inv(G)
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)
Irinj = (Pinj*V[buspi-1].real + Qinj*V[busqi-1].imag)/(numpy.absolute(V[buspi-1])**2)
Ixinj = (Pinj*V[buspi-1].imag - Qinj*V[busqi-1].real)/(numpy.absolute(V[buspi-1])**2)
z[pidx] = Irinj
z[qidx] = Ixinj
Irbr = (Pbr*V.real[fbuspf-1] + Qbr*V.imag[fbusqf-1])/(V.mag[fbuspf-1]**2)
Ixbr = (Pbr*V.imag[fbuspf-1] - Qbr*V.real[fbusqf-1])/(V.mag[fbuspf-1]**2)
Irbr = (Pbr*V[fbuspf-1].real + Qbr*V[fbusqf-1].imag)/(numpy.absolute(V[fbuspf-1])**2)
Ixbr = (Pbr*V[fbuspf-1].imag - Qbr*V[fbusqf-1].real)/(numpy.absolute(V[fbuspf-1])**2)
z[pfidx] = Irbr
z[qfidx] = Ixbr
......@@ -487,15 +488,12 @@ def DssePmu(system, zdata, Ymatrix, Adj):
V.real = State[:nodes_num]
V.imag = State[nodes_num:]
V = Real_to_all(V.real, V.imag)
num_iter = num_iter+1
return V
def DsseMixed(system, zdata, Ymatrix, Adj):
def DsseMixed(system, zdata):
""" It performs state estimation using rectangular node voltage state variables
and it is built to work in scenarios where both conventional and PMU measurements
are simultaneously present."""
......@@ -544,11 +542,10 @@ def DsseMixed(system, zdata, Ymatrix, Adj):
weights = zdata.mstddev**(-2)
W = numpy.diag(weights)
Admittance = Complex_to_all(Ymatrix)
Gmatrix = Admittance.real
Bmatrix = Admittance.imag
Yabs_matrix = Admittance.mag
Yphase_matrix = Admittance.phase
Gmatrix = system.Ymatrix.real
Bmatrix = system.Ymatrix.imag
Yabs_matrix = numpy.absolute(system.Ymatrix)
Yphase_matrix = numpy.angle(system.Ymatrix)
""" Jacobian for Power Injection Measurements (converted to equivalent
rectangualar current measurements) """
......@@ -561,7 +558,7 @@ def DsseMixed(system, zdata, Ymatrix, Adj):
H2[i][m2] = Bmatrix[m][m]
H3[i][m] = - Bmatrix[m][m]
H3[i][m2] = - Gmatrix[m][m]
idx = numpy.subtract(Adj[m],1)
idx = numpy.subtract(system.Adj[m],1)
H2[i][idx] = - Gmatrix[m][idx]
H3[i][idx] = - Bmatrix[m][idx]
idx2 = idx + len(system.nodes)
......@@ -641,7 +638,8 @@ def DsseMixed(system, zdata, Ymatrix, Adj):
epsilon = 5
Vr = numpy.ones(len(system.nodes))
Vx = numpy.zeros(len(system.nodes))
V = Real_to_all(Vr, Vx)
V = Vr + 1j*Vx
#V = Real_to_all(Vr, Vx)
num_iter = 0
StateVr = numpy.ones(len(system.nodes))
......@@ -650,24 +648,24 @@ def DsseMixed(system, zdata, Ymatrix, Adj):
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)
Irinj = (Pinj*V[buspi-1].real + Qinj*V[busqi-1].imag)/(numpy.absolute(V[buspi-1])**2)
Ixinj = (Pinj*V[buspi-1].imag - Qinj*V[busqi-1].real)/(numpy.absolute(V[buspi-1])**2)
z[pidx] = Irinj
z[qidx] = Ixinj
Irbr = (Pbr*V.real[fbuspf-1] + Qbr*V.