Commit 0a0265c5 authored by martin.moraga's avatar martin.moraga
Browse files

fixed errors

parent 8f8c7add
......@@ -24,7 +24,7 @@ class Measurement():
def __init__(self, element, element_type, meas_type, meas_value, unc):
"""
Creates a measurement, which is used by the estimation module. Possible types of measurements are: v, p, q, i, Vpmu and Ipmu
@element: pointer to measured element
@element: pointer to the topology_node / topology_branch (object of class network.Node / network.Branch)
@element_type: Clarifies which element is measured.
@meas_type:
@meas_value: measurement value.
......@@ -42,12 +42,6 @@ class Measurement():
self.meas_type = meas_type
self.meas_value = meas_value
self.std_dev = unc/300
"""
if meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
self.std_dev = meas_value*unc/100
elif meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
self.std_dev = unc/100
"""
self.mval = 0.0 #measured values (affected by uncertainty)
class Measurents_set():
......@@ -64,7 +58,7 @@ class Measurents_set():
"""
read measurements from file.
@param powerflow_results
@param powerflow_results: object of class acs.state_estimation.results.Results
@param file_name
"""
with open(file_name) as json_file:
......@@ -72,97 +66,111 @@ class Measurents_set():
for key, value in data['Measurement'].items():
if key=="Vmag":
unc = float(value['unc'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_node
meas_value = np.abs(powerflow_results.nodes[index-1].voltage)
for uuid in value['uuid']:
pf_node = powerflow_results.get_node(uuid=uuid)
element = pf_node.topology_node
meas_value = np.abs(pf_node.voltage_pu)
self.create_measurement(element, ElemType.Node, MeasType.V_mag, meas_value, unc)
elif key=="Imag":
unc = float(value['unc'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_branch
meas_value =np.abs(powerflow_results.branches[index-1].current)
for uuid in value['uuid']:
pf_branch = powerflow_results.get_branch(uuid=uuid)
element = pf_branch.topology_branch
meas_value =np.abs(pf_branch.current_pu)
self.create_measurement(element, ElemType.Branch, MeasType.I_mag, meas_value, unc)
elif key=="Pinj":
unc = float(value['unc'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_node
meas_value = powerflow_results.nodes[index-1].power.real
for uuid in value['uuid']:
pf_node = powerflow_results.get_node(uuid=uuid)
element = pf_node.topology_node
meas_value = pf_node.power_pu.real
self.create_measurement(element, ElemType.Node, MeasType.Sinj_real, meas_value, unc)
elif key=="Qinj":
unc = float(value['unc'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_node
meas_value = powerflow_results.nodes[index-1].power.imag
for uuid in value['uuid']:
pf_node = powerflow_results.get_node(uuid=uuid)
element = pf_node.topology_node
meas_value = pf_node.power_pu.imag
self.create_measurement(element, ElemType.Node, MeasType.Sinj_imag, meas_value, unc)
elif key=="P1":
unc = float(value['unc'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_branch
meas_value = powerflow_results.branches[index-1].power.real
for uuid in value['uuid']:
pf_branch = powerflow_results.get_branch(uuid=uuid)
element = pf_branch.topology_branch
meas_value = pf_branch.power_pu.real
self.create_measurement(element, ElemType.Branch, MeasType.S1_real, meas_value, unc)
elif key=="Q1":
unc = float(value['unc'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_branch
meas_value = powerflow_results.branches[index-1].power.imag
for uuid in value['uuid']:
pf_branch = powerflow_results.get_branch(uuid=uuid)
element = pf_branch.topology_branch
meas_value = pf_branch.power_pu.imag
self.create_measurement(element, ElemType.Branch, MeasType.S1_imag, meas_value, unc)
elif key=="P2":
unc = float(value['unc'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_branch
meas_value = powerflow_results.branches[index-1].power2.real
for uuid in value['uuid']:
pf_branch = powerflow_results.get_branch(uuid=uuid)
element = pf_branch.topology_branch
meas_value = pf_branch.power2_pu.real
self.create_measurement(element, ElemType.Branch, MeasType.S2_real, meas_value, unc)
elif key=="Q2":
unc = float(value['unc'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_branch
meas_value = powerflow_results.branches[index-1].power2.imag
for uuid in value['uuid']:
pf_branch = powerflow_results.get_branch(uuid=uuid)
element = pf_branch.topology_branch
meas_value = pf_branch.power2_pu.imag
self.create_measurement(element, ElemType.Branch, MeasType.S2_imag, meas_value, unc)
elif key=="Vpmu":
unc_mag = float(value['unc_mag'])
unc_phase = float(value['unc_phase'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_node
meas_value_mag = np.abs(powerflow_results.nodes[index-1].voltage)
#meas_value_phase = np.angle(powerflow_results.nodes[index-1].voltage)
for uuid in value['uuid']:
pf_node = powerflow_results.get_node(uuid=uuid)
element = pf_node.topology_node
meas_value_mag = np.abs(pf_node.voltage_pu)
self.create_measurement(element, ElemType.Node, MeasType.Vpmu_mag, meas_value_mag, unc_mag)
#self.create_measurement(element, ElemType.Node, MeasType.Vpmu_phase, meas_value_phase, unc_phase)
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_node
meas_value_phase = np.angle(powerflow_results.nodes[index-1].voltage)
for uuid in value['uuid']:
pf_node = powerflow_results.get_node(uuid=uuid)
element = pf_node.topology_node
meas_value_phase = np.angle(pf_node.voltage_pu)
self.create_measurement(element, ElemType.Node, MeasType.Vpmu_phase, meas_value_phase, unc_phase)
elif key=="Ipmu":
unc_mag = float(value['unc_mag'])
unc_phase = float(value['unc_phase'])
for index in value['idx']:
element = powerflow_results.nodes[index-1].topology_branch
meas_value_mag =np.abs(powerflow_results.branches[index-1].current)
meas_value_phase =np.angle(powerflow_results.branches[index-1].current)
for uuid in value['uuid']:
pf_branch = powerflow_results.get_branch(uuid=uuid)
element = pf_branch.topology_branch
meas_value_mag = np.abs(pf_branch.current_pu)
meas_value_phase = np.angle(pf_branch.current_pu)
self.create_measurement(element, ElemType.Branch, MeasType.Ipmu_mag, meas_value_mag, unc_mag)
self.create_measurement(element, ElemType.Branch, MeasType.Ipmu_phase, meas_value_phase, unc_phase)
def meas_creation(self, seed=None):
def meas_creation(self, dist="normal", seed=None):
"""
It calculates the measured values (affected by uncertainty) at the measurement points
which distribution should be used? if gaussian --> stddev must be divided by 3
@param seed: Seed for RandomState (to make the random numbers predictable)
Must be convertible to 32 bit unsigned integers.
@param seed: - normal: use normal distribution (-->std_dev are divided by 3)
- uniform: use normal distribution
"""
if seed is None:
np.random.seed(seed)
if dist == "normal":
err_pu = np.random.normal(0,1,len(self.measurements))
#err_pu = np.random.uniform(-1,1,len(self.measurements))
else:
np.random.seed(seed)
for index, measurement in enumerate(self.measurements):
if measurement.meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.meas_value*measurement.std_dev
elif measurement.meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.std_dev
measurement.mval = measurement.meas_value + zdev*err_pu[index]
elif dist == "uniform":
err_pu = np.random.uniform(-1,1,len(self.measurements))
for index, measurement in enumerate(self.measurements):
if measurement.meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.meas_value*measurement.std_dev
elif measurement.meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.std_dev
measurement.mval = measurement.meas_value + zdev*err_pu[index]
for index, measurement in enumerate(self.measurements):
if measurement.meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = (measurement.meas_value*measurement.std_dev)
elif measurement.meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.std_dev
measurement.mval = measurement.meas_value + np.multiply(3*zdev,err_pu[index])
def meas_creation_test(self, err_pu):
"""
......@@ -171,7 +179,12 @@ class Measurents_set():
This function takes as paramenter a random distribution.
"""
for index, measurement in enumerate(self.measurements):
if measurement.meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
measurement.std_dev = np.absolute(measurement.meas_value)*measurement.std_dev
elif measurement.meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
measurement.std_dev = measurement.std_dev
measurement.mval = measurement.meas_value + measurement.std_dev*err_pu[index]
#measurement.mval = measurement.meas_value + measurement.std_dev*err_pu[index]
def getMeasurements(self, type):
"""
......@@ -217,22 +230,10 @@ class Measurents_set():
#the weight is small and can bring instability during matrix inversion, so we "cut" everything below 10^-6
if measurement.std_dev<10**(-6):
measurement.std_dev = 10**(-6)
#weights[index] = (measurement.std_dev/3)**(-2)
weights[index] = (measurement.std_dev)**(-2)
return weights
def getMeasValues(self):
"""
for test purposes
returns an array with all measured values
"""
meas_val = np.zeros(len(self.measurements))
for index, measurement in enumerate(self.measurements):
meas_val[index] = measurement.meas_value
return meas_val
def getmVal(self):
"""
returns an array with all measured values (affected by uncertainty)
......@@ -276,12 +277,39 @@ class Measurents_set():
return std_dev
def getmVal_test(self):
def getMeasValues(self, type=None):
"""
for test purposes
returns an array with all measured values
"""
if type is None:
meas_val = np.zeros(len(self.measurements))
for index, measurement in enumerate(self.measurements):
meas_val[index] = measurement.meas_value
else:
meas_val = np.zeros(self.getNumberOfMeasurements(type=type))
idx=0
for index, measurement in enumerate(self.measurements):
if measurement.meas_type is type:
meas_val[idx] = measurement.meas_value
idx += 1
return meas_val
def getmVal_test(self, type=None):
"""
returns an array with all measured values (affected by uncertainty)
"""
mVal = np.zeros(len(self.measurements))
for index, measurement in enumerate(self.measurements):
mVal[index] = measurement.mval
if type is None:
mVal = np.zeros(len(self.measurements))
for index, measurement in enumerate(self.measurements):
mVal[index] = measurement.mval
else:
mVal = np.zeros(self.getNumberOfMeasurements(type=type))
idx=0
for index, measurement in enumerate(self.measurements):
if measurement.meas_type is type:
mVal[idx] = measurement.mval
idx += 1
return mVal
\ No newline at end of file
......@@ -10,25 +10,27 @@ class BusType(Enum):
pq = 3
class Node():
def __init__(self, name='', uuid='', base_voltage = 1.0, base_apparent_power=1.0,
v_mag=0.0, v_phase=0.0, p=0.0, q=0.0, index=0, bus_type='PQ'):
self.name = name
def __init__(self, uuid='', base_voltage = 1.0, base_apparent_power=1.0, v_mag=0.0,
v_phase=0.0, p=0.0, q=0.0, index=0, bus_type='PQ', pu=False):
self.uuid = uuid
self.index = index
self.baseVoltage = base_voltage
self.base_apparent_power = base_apparent_power
self.base_current = self.base_apparent_power/self.baseVoltage
self.type = BusType[bus_type]
self.voltage = v_mag*np.cos(v_phase) + 1j * v_mag*np.sin(v_phase)
self.power = complex(p, q)
self.power_pu = complex(p, q)/self.base_apparent_power
self.voltage_pu = self.voltage/self.baseVoltage
class Branch():
def __init__(self, r, x, start_node, end_node, base_voltage=1.0, base_apparent_power=1.0):
def __init__(self, uuid='', r=0.0, x=0.0, start_node=None, end_node=None,
base_voltage=1.0, base_apparent_power=1.0):
self.uuid = uuid
self.baseVoltage = base_voltage
self.base_apparent_power = base_apparent_power
self.base_impedance = base_voltage**2 / self.base_apparent_power
self.base_current = self.base_apparent_power/self.baseVoltage
self.base_impedance = base_voltage**2/self.base_apparent_power
self.start_node = start_node
self.end_node = end_node
self.r = r
......@@ -44,10 +46,6 @@ class System():
def __init__(self):
self.nodes=[]
self.branches=[]
self.bR=[]
self.bX=[]
self.P=[]
self.Q=[]
self.Ymatrix = np.zeros([], dtype=np.complex)
self.Adjacencies = np.array([])
......@@ -55,7 +53,7 @@ class System():
"""
To fill the vectors node, branch, bR, bX, P and Q
"""
index = 0
for uuid, element in res.items():
if element.__class__.__name__=="TopologicalNode":
vmag = 0.0
......@@ -73,19 +71,15 @@ class System():
if element2.getNodeUUID() == uuid:
pInj += element2.p
qinj += element2.q
self.P.append(pInj)
self.Q.append(qinj)
index=len(self.P)-1
node_type = self._getNodeType(element)
base_voltage = element.BaseVoltage.nominalVoltage
self.nodes.append(Node(name=element.name, uuid=uuid, base_voltage=base_voltage,
base_apparent_power=base_apparent_power, v_mag=vmag,
v_phase=vphase, p=pInj, q=qinj, index=index, bus_type=node_type))
self.nodes.append(Node(uuid=uuid, base_voltage=base_voltage, v_mag=vmag,
base_apparent_power=base_apparent_power, v_phase=vphase,
p=pInj, q=qinj, index=index, bus_type=node_type))
index = index+1
for uuid, element in res.items():
if element.__class__.__name__=="ACLineSegment":
self.bR.append(element.r)
self.bX.append(element.x)
for node in self.nodes:
if element.startNodeID==node.uuid:
startNode = node
......@@ -94,17 +88,12 @@ class System():
if element.endNodeID==node.uuid:
endNode=node
break
base_voltage = element.BaseVoltage.nominalVoltage
self.branches.append(Branch(r=element.r, x=element.x, start_node=startNode,
self.branches.append(Branch(uuid=uuid, r=element.r, x=element.x, start_node=startNode,
end_node=endNode, base_voltage=base_voltage,
base_apparent_power=base_apparent_power))
continue
elif element.__class__.__name__=="PowerTransformer":
bR = element.primaryConnection.r
bX = element.primaryConnection.x
self.bR.append(bR)
self.bX.append(bX)
for i in range(len(self.nodes)):
if element.startNodeID==self.nodes[i].uuid:
startNode=self.nodes[i]
......@@ -113,16 +102,16 @@ class System():
if element.endNodeID==self.nodes[i].uuid:
endNode=self.nodes[i]
break
#base voltage = high voltage side (=primaryConnection)
base_voltage = element.primaryConnection.BaseVoltage.nominalVoltage
self.branches.append(Branch(r=bR, x=bX, start_node=startNode, end_node=endNode,
base_voltage=base_voltage, base_apparent_power=base_apparent_power))
self.branches.append(Branch(uuid=uuid, r=element.primaryConnection.r, x=element.primaryConnection.x,
start_node=startNode, end_node=endNode, base_voltage=base_voltage,
base_apparent_power=base_apparent_power))
continue
else:
continue
#calculate impedance matrix
#calculate admitance matrix
self.Ymatrix_calc()
def _getNodeType(self, node):
......@@ -150,23 +139,4 @@ class System():
self.Ymatrix[fr][fr] += branch.y_pu
self.Ymatrix[to][to] += branch.y_pu
self.Adjacencies[fr].append(to+1) #to + 1???
self.Adjacencies[to].append(fr+1) #fr + 1???
def load_python_data(nodes, branches, type):
system = System()
for node_idx in range(0, nodes.num):
if BusType[type[node_idx]] == BusType.slack:
system.nodes.append(Node(v_mag=nodes.P2[0], v_phase=nodes.Q2[0], p=0, q=0, index=node_idx, bus_type=type[node_idx]))
elif BusType[type[node_idx]] == BusType.PQ:
system.nodes.append(Node(v_mag=0, v_phase=0, p=nodes.P2[node_idx], q=nodes.Q2[node_idx], index=node_idx, bus_type=type[node_idx]))
elif BusType[type[node_idx]] == BusType.PV:
#TODO
pass
for branch_idx in range(0, branches.num):
system.branches.append(Branch(branches.R[branch_idx], branches.X[branch_idx],
system.nodes[branches.start[branch_idx]-1], system.nodes[branches.end[branch_idx]-1]))
system.Ymatrix_calc()
return system
\ No newline at end of file
self.Adjacencies[to].append(fr+1) #fr + 1???
\ No newline at end of file
import sys
import math
import numpy as np
from .network import BusType
from .results import Results
......@@ -10,20 +8,20 @@ def solve(system):
nodes_num = len(system.nodes)
branches_num = len(system.branches)
z = np.zeros(2*(nodes_num))
h = np.zeros(2*(nodes_num))
H = np.zeros((2*(nodes_num),2*(nodes_num)))
z = np.zeros(2*nodes_num)
h = np.zeros(2*nodes_num)
H = np.zeros((2*nodes_num,2*nodes_num))
for i in range(0,nodes_num):
m = 2*i
i2 = i + nodes_num
type = system.nodes[i].type
if type is BusType.SLACK:
node_type = system.nodes[i].type
if node_type == BusType.SLACK:
z[m] = np.real(system.nodes[i].voltage_pu)
z[m+1] = np.imag(system.nodes[i].voltage_pu)
H[m][i] = 1
H[m+1][i2] = 1
elif type is BusType.PQ:
elif node_type is BusType.PQ:
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])
......@@ -34,7 +32,7 @@ def solve(system):
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:
elif node_type is BusType.PV:
z[m+1] = np.real(system.nodes[i].power)
H[m][i] = - np.real(system.Ymatrix[i][i])
H[m][i2] = np.imag(system.Ymatrix[i][i])
......@@ -43,8 +41,7 @@ def solve(system):
H[m][idx1] = - np.real(system.Ymatrix[i][idx1])
H[m][idx2] = np.imag(system.Ymatrix[i][idx1])
#epsilon = 10**(-10)
epsilon = 10**(-3)
epsilon = 10**(-10)
diff = 5
V = np.ones(nodes_num) + 1j* np.zeros(nodes_num)
num_iter = 0
......@@ -56,22 +53,22 @@ def solve(system):
for i in range(0, nodes_num):
m = 2*i
i2 = i + nodes_num
type = system.nodes[i].type
if type is BusType.SLACK:
node_type = system.nodes[i].type
if node_type is BusType.SLACK:
h[m] = np.inner(H[m],State)
h[m+1] = np.inner(H[m+1],State)
elif type is BusType.PQ:
elif node_type is BusType.PQ:
z[m] = (np.real(system.nodes[i].power_pu)*np.real(V[i]) + np.imag(system.nodes[i].power_pu)*np.imag(V[i]))/(np.abs(V[i])**2)
z[m+1] = (np.real(system.nodes[i].power_pu)*np.imag(V[i]) - np.imag(system.nodes[i].power_pu)*np.real(V[i]))/(np.abs(V[i])**2)
h[m] = np.inner(H[m],State)
h[m+1] = np.inner(H[m+1],State)
elif type is BusType.PV:
elif node_type is BusType.PV:
z[m] = (np.real(system.nodes[i].power_pu)*np.real(V[i])+ np.imag(system.nodes[i].power_pu)*np.imag(V[i]))(np.abs(V[i])**2)
h[m] = np.inner(H[m],State)
h[m+1] = np.abs(V[i])
H[m+1][i] = np.cos(np.angle(V[i]))
H[m+1][i2] = np.sin(np.angle(V[i]))
r = np.subtract(z,h)
Hinv = np.linalg.inv(H)
Delta_State = np.inner(Hinv,r)
......
import numpy as np
from results import Results
from measurement import *
from acs.state_estimation.results import Results
from acs.state_estimation.measurement import *
def DsseCall(system, measurements):
"""
......@@ -211,7 +211,6 @@ def DssePmu(nodes_num, measurements, Gmatrix, Bmatrix, Adj):
# Iteration of Netwon Rapson method: needed to solve non-linear system of equation
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
z = convertSinjMeasIntoCurrents(measurements, V, z, pidx, qidx)
......@@ -230,7 +229,6 @@ def DssePmu(nodes_num, measurements, Gmatrix, Bmatrix, Adj):
V.real = State[:nodes_num]
V.imag = State[nodes_num:]
#print(V.real)
num_iter = num_iter+1
return V
......@@ -296,14 +294,12 @@ def DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_mat
epsilon = 5
num_iter = 0
k = 0
# Iteration of Netwon Rapson method: needed to solve non-linear system of equation
while epsilon>10**(-6):
k = k+1
""" Computation of equivalent current measurements in place of the power measurements """
z = convertSinjMeasIntoCurrents(measurements, V, z, pidx, qidx)
z = convertSbranchMeasIntoCurrents(measurements, V, z, p1br, q1br, p2br, q2br)
""" Voltage Magnitude Measurements """
h1, H1 = update_h1_vector(measurements, V, vidx, nvi, nodes_num, type=1)
......@@ -339,15 +335,12 @@ def DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_mat
State = State + Delta_State
epsilon = np.amax(np.absolute(Delta_State))
V.real = State[:nodes_num]
V.imag = State[nodes_num:]
num_iter = num_iter+1
#if k==4:
# return H
return V
def calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, Adj, type):
......@@ -612,7 +605,7 @@ def update_h1_vector(measurements, V, vidx, nvi, nodes_num, type):
H1[i][m2] = np.sin(np.angle(V[node_index]))
elif type==2:
if m > 0:
m2 = node_index + node.num - 1
m2 = node_index + nodes_num - 1
H1[i][m2] = np.sin(V.phase[m])
return h1, H1
......
import numpy as np
import cmath
import sys
import numpy as np
from villas.dataprocessing.readtools import read_timeseries_dpsim
class ResultsNode():
......@@ -9,14 +8,20 @@ class ResultsNode():
self.voltage = complex(0, 0)
self.current = complex(0, 0)
self.power = complex(0, 0)
self.voltage_pu = complex(0, 0)
self.current_pu = complex(0, 0)
self.power_pu = complex(0, 0)
class ResultsBranch():
def __init__(self, topo_branch):
self.topology_branch = topo_branch
self.current = 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