fixed some errors

acs/state_estimation/measurement.py

 from enum import Enum
+import json
 import numpy as np
 
 class ElemType(Enum):
@@ -20,14 +21,14 @@ class MeasType(Enum):
 	S2_imag = 12						#Reactive Power flow at branch, measured at final node (S2.imag)
 
 class Measurement():
-	def __init__(self, element, element_type, meas_type, meas_value, std_dev):
+	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_type: Clarifies which element is measured.
 		@meas_type: 
 		@meas_value: measurement value.
-		@std_dev: standard deviation in the same unit as the measurement.
+		@unc: Measurement uncertainty in percent
 		"""
 		
 		if not isinstance(element_type, ElemType):
@@ -40,33 +41,123 @@ class Measurement():
 		self.element_type = element_type
 		self.meas_type = meas_type
 		self.meas_value = meas_value
-		self.std_dev = std_dev
-		self.std_dev = std_dev
+		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():
 	def __init__(self):
 		self.measurements = []			#array with all measurements
 		
-	def create_measurement(self, element, element_type, meas_type, meas_value, std_dev):
+	def create_measurement(self, element, element_type, meas_type, meas_value, unc):
 		"""
 		to add elements to the measurements array
 		"""
-		self.measurements.append(Measurement(element, element_type, meas_type, meas_value, std_dev))
+		self.measurements.append(Measurement(element, element_type, meas_type, meas_value, unc))
 	
-	def meas_creation(self):
+	def read_measurements_from_file(self, powerflow_results, file_name):
+		"""
+		read measurements from file.
+
+		@param powerflow_results 
+		@param file_name
+		"""
+		with open(file_name) as json_file:
+			data = json.load(json_file)
+
+		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)
+					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)
+					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
+					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
+					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
+					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
+					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
+					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
+					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)
+					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)
+			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)
+					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):
 		""" 
 		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.
 		"""
-		err_pu = np.random.normal(0,1,len(self.measurements))
+		if seed is None:
+			#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)
+			err_pu = np.random.uniform(-1,1,len(self.measurements))
+		
 		for index, measurement in enumerate(self.measurements):
-			measurement.mval = measurement.meas_value + self.std_dev*err_pu[index]
+			measurement.mval = measurement.meas_value + measurement.std_dev*err_pu[index]
 
 	def meas_creation_test(self, err_pu):
 		""" 
 		For test purposes.
 		It calculates the measured values (affected by uncertainty) at the measurement points.
-		This function takes as paramenter the random gaussian distribution. 
+		This function takes as paramenter a random distribution. 
 		"""
 		for index, measurement in enumerate(self.measurements):
 			measurement.mval = measurement.meas_value + measurement.std_dev*err_pu[index]

acs/state_estimation/network.py

@@ -39,6 +39,7 @@ class System():
 		self.Ymatrix = np.zeros((1, 1),dtype=np.complex)
 		self.Adjacencies = np.array([])
 
+	#def load_cim_data(self, res, Sb, Vb, Zb):
 	def load_cim_data(self, res):
 		#this function is used to fill the vectors node, branch, bR, bX, P and Q
 		for key, value in res.items():
@@ -77,11 +78,17 @@ class System():
 				self.branches.append(Branch(value.primaryConnection.r, value.primaryConnection.x, startNode, endNode))
 			else:
 				continue
+
+		#determine the impedance matrix
+		#self.Ymatrix_calc(Zb)
 		self.Ymatrix_calc()
 
+	#def Ymatrix_calc(self, Zb):
 	def Ymatrix_calc(self):
+		"""
+		@param Zb: base value of impedance
+		"""
 		nodes_num = len(self.nodes)
-		branches_num = len(self.branches)
 		self.Ymatrix = np.zeros((nodes_num, nodes_num),dtype=np.complex)
 		self.Adjacencies = [[] for _ in range(nodes_num)]
 		for branch in self.branches:
@@ -94,6 +101,8 @@ class System():
 			self.Adjacencies[fr].append(to+1)	#to + 1???
 			self.Adjacencies[to].append(fr+1)	#fr + 1???
 
+		#self.Ymatrix = self.Ymatrix*Zb
+
 def load_python_data(nodes, branches, type):
 	system = System()
 	
@@ -108,20 +117,4 @@ def load_python_data(nodes, branches, type):
 		system.nodes[branches.start[branch_idx]-1], system.nodes[branches.end[branch_idx]-1]))
 	
 	system.Ymatrix_calc()
-	return system
-	
-def Ymatrix_calc(system):
-	nodes_num = len(system.nodes)
-	branches_num = len(system.branches)
-	Ymatrix = np.zeros((nodes_num, nodes_num),dtype=np.complex)
-	Adjacencies = [[] for _ in range(nodes_num)]
-	for index in range(branches_num):
-		fr = system.branches[index].start_node.index
-		to = system.branches[index].end_node.index
-		Ymatrix[fr][to] -= system.branches[index].y
-		Ymatrix[to][fr] -= system.branches[index].y
-		Ymatrix[fr][fr] += system.branches[index].y
-		Ymatrix[to][to] += system.branches[index].y
-		Adjacencies[fr].append(to+1)
-		Adjacencies[to].append(fr+1)
-	return Ymatrix, Adjacencies
+	return system
\ No newline at end of file

acs/state_estimation/nv_powerflow_cim.py

 import sys
 import math
 import numpy as np
-from network import Ymatrix_calc
 from network import BusType
 
 sys.path.append("../../../dataprocessing")

acs/state_estimation/nv_state_estimator_cim.py

-import sys
 import numpy as np
 from results import Results
 from measurement import *
@@ -41,19 +40,14 @@ def DsseCall(system, measurements):
 	if est_code == 1:
 		Vest = DsseTrad(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix)
 	elif est_code == 2:
-		Vest = DssePmu(nodes_num, measurements, Gmatrix, Bmatrix)
+		Vest = DssePmu(nodes_num, measurements, Gmatrix, Bmatrix, Adj)
 	else:
 		Vest = DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix, Adj)
 
 	# calculate all the other quantities of the grid
 	results = Results(system)
 	results.load_voltages(Vest)
-	results.calculateI()
-	results.calculateIinj()
-	results.calculateSinj()
-	results.calculateI()
-	results.calculateS1()
-	results.calculateS2()
+	results.calculate_all()
 
 	return results
 
@@ -114,7 +108,7 @@ def DsseTrad(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matr
 		z = convertSbranchMeasIntoCurrents(measurements, V, z, p1br, q1br, p2br, q2br)
 
 		""" Voltage Magnitude Measurements """
-		h1, H1 = update_h1_vector(measurements, V, vidx, nvi, nodes_num, type==2)
+		h1, H1 = update_h1_vector(measurements, V, vidx, nvi, nodes_num, type=2)
 				
 		""" Power Injection Measurements """
 		# h(x) vector where power injections are present
@@ -127,7 +121,7 @@ def DsseTrad(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matr
 		h5 = np.inner(H5,State)
 		
 		""" Current Magnitude Measurements """
-		h6, H6 = update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nodes_num, num_iter, type==2)
+		h6, H6 = update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nodes_num, num_iter, type=2)
 		
 		""" WLS computation """
 		# all the sub-matrixes of H calcualted so far are merged in a unique matrix
@@ -177,7 +171,7 @@ def DssePmu(nodes_num, measurements, Gmatrix, Bmatrix, Adj):
 	W = np.diag(weights)
 	
 	#Jacobian for Power Injection Measurements
-	H2, H3 = calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, Adj, type==1)
+	H2, H3 = calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, Adj, type=1)
 
 	#Jacobian for branch Power Measurements
 	H4, H5 = calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type=1)
@@ -265,7 +259,7 @@ def DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_mat
 	# Jacobian Matrix. Includes the derivatives of the measurements (voltages, currents, powers) with respect to the states (voltages)
 
 	#Jacobian for Power Injection Measurements
-	H2, H3 = calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, Adj, type==1)
+	H2, H3 = calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, Adj, type=1)
 
 	#Jacobian for branch Power Measurements
 	H4, H5 = calculateJacobiBranchPower(measurements, nodes_num, Gmatrix, Bmatrix, type=1)
@@ -311,7 +305,7 @@ def DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_mat
 		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)
+		h1, H1 = update_h1_vector(measurements, V, vidx, nvi, nodes_num, type=1)
 
 		""" Power Injection Measurements """
 		h2 = np.inner(H2,State)
@@ -322,7 +316,7 @@ def DsseMixed(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_mat
 		h5 = np.inner(H5,State)
 
 		""" Current Magnitude Measurements """
-		h6, H6 = update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nodes_num, num_iter, type==1)
+		h6, H6 = update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nodes_num, num_iter, type=1)
 		
 		""" PMU Voltage Measurements """
 		h7 = np.inner(H7,State)
@@ -374,7 +368,6 @@ def calculateJacobiMatrixSinj(measurements, nodes_num, Gmatrix, Bmatrix, Adj, ty
 	pinj_meas = measurements.getMeasurements(type=MeasType.Sinj_real)
 	#get all measurements of type MeasType.Sinj_real
 	qinj_meas = measurements.getMeasurements(type=MeasType.Sinj_imag)	
-	
 	if type == 1:
 		H2 = np.zeros((len(pinj_meas), 2*nodes_num))
 		H3 = np.zeros((len(qinj_meas), 2*nodes_num))
@@ -504,6 +497,7 @@ def calculateJacobiVoltagePmu(measurements, nodes_num, Gmatrix, Bmatrix):
 	H7 = np.zeros((len(Vpmu_mag_meas),2*nodes_num))
 	H8 = np.zeros((len(Vpmu_mag_meas),2*nodes_num))
 	
+	#TODO: index of Vmag = index of Vphase???
 	for index, measurement in enumerate(Vpmu_mag_meas):
 		vamp = measurement.meas_value
 		vtheta = Vpmu_phase_meas[index].meas_value
@@ -668,11 +662,11 @@ def update_h6_vector(measurements, V, iidx, nii, Yabs_matrix, Yphase_matrix, nod
 			H6[i][n2] = Yabs_matrix[m][n]*(np.cos(Yphase_matrix[m][n])*h6im[i] - np.sin(Yphase_matrix[m][n])*h6re[i])/h6[i]
 		if type==2:
 			if m > 0:
-                m2 = m + node.num-1
-                H6[i][m2] = - Yabs_matrix[m][n]*(np.cos(Yphase_matrix[m][n])*h6im[i] - np.sin(Yphase_matrix[m][n])*h6re[i])/h6[i]
-            if n > 0:
-                n2 = n + node.num-1
-                H6[i][n2] = Yabs_matrix[m][n]*(np.cos(Yphase_matrix[m][n])*h6im[i] - np.sin(Yphase_matrix[m][n])*h6re[i])/h6[i]
+				m2 = m + node.num-1
+				H6[i][m2] = - Yabs_matrix[m][n]*(np.cos(Yphase_matrix[m][n])*h6im[i] - np.sin(Yphase_matrix[m][n])*h6re[i])/h6[i]
+			if n > 0:
+				n2 = n + node.num-1
+				H6[i][n2] = Yabs_matrix[m][n]*(np.cos(Yphase_matrix[m][n])*h6im[i] - np.sin(Yphase_matrix[m][n])*h6re[i])/h6[i]
 			
 	return h6, H6
 

acs/state_estimation/results.py

 import numpy as np
+import cmath
 
 class ResultsNode():
 	def __init__(self, topo_node):		
@@ -41,7 +42,27 @@ class Results():
 			for elem in self.nodes:
 				if elem.topology_node.index == index:
 					elem.voltage = V[index]
-					
+
+	def load_voltages_from_dict(self, V):
+		"""
+		load the voltages of V-dict 
+		"""
+		for key, value in V.items():
+			for elem in self.nodes:
+				if elem.topology_node.index == index:
+					elem.voltage = V[index]
+
+	def calculate_all(self):
+		"""
+		calculate all quantities of the grid
+		"""
+		self.calculateI()
+		self.calculateIinj()
+		self.calculateSinj()
+		self.calculateIinj()
+		self.calculateS1()
+		self.calculateS2()
+
 	def calculateI(self):
 		"""
 		To calculate the branch currents
@@ -92,13 +113,21 @@ class Results():
 				if branch_index == node.topology_node.index:
 					branch.power2 = -node.voltage*(np.conj(branch.current))
 	
-	def get_node(self, index):
+	def get_node(self, index=None, uuid=None):
 		"""
 		return the PowerflowNode with PowerflowNode.topology_node.index == index
+		returns a node with a certain uuid or a certain index (not both!):
+		- if index in not None --> return the PowerflowNode with PowerflowNode.topology_node.index == index
+		- if uuid in not None --> return the PowerflowNode with PowerflowNode.topology_node.uuid == uuid
 		"""
-		for node in self.nodes:
-			if index == node.topology_node.index:
-				return node
+		if index is not None:
+			for node in self.nodes:
+				if index == node.topology_node.index:
+					return node
+		elif uuid is not None:
+			for node in self.nodes:
+				if uuid == node.topology_node.uuid:
+					return node
 				
 	def get_voltages(self):
 		"""
@@ -159,4 +188,11 @@ class Results():
 		S2 = np.array([])
 		for branch in self.branches:
 			S2 = np.append(S2, branch.power2)
-		return S2
\ No newline at end of file
+		return S2
+
+	def print_voltages_polar(self):
+		"""
+		for test purposes
+		"""
+		for node in self.nodes:
+			print(node.topology_node.uuid + ": " + str(cmath.polar(node.voltage)))
\ No newline at end of file