improve network.py and fixed some errors

acs/state_estimation/measurement.py

@@ -41,7 +41,7 @@ class Measurement():
 		self.element_type = element_type
 		self.meas_type = meas_type
 		self.meas_value = meas_value
-		self.std_dev = unc/100
+		self.std_dev = unc/300
 		"""
 		if meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
 			self.std_dev = meas_value*unc/100
@@ -150,8 +150,8 @@ class Measurents_set():
 					 Must be convertible to 32 bit unsigned integers.
 		"""
 		if seed is None:
-			#err_pu = np.random.normal(0,1,len(self.measurements))
-			err_pu = np.random.uniform(-1,1,len(self.measurements))
+			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))
@@ -217,8 +217,9 @@ 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/3)**(-2)
+			weights[index] = (measurement.std_dev)**(-2)
+
 		return weights
 	
 	def getMeasValues(self):

acs/state_estimation/network.py

@@ -9,24 +9,36 @@ class BusType(Enum):
 	PQ = 3
 	pq = 3
 
-
 class Node():
-	def __init__(self, name='', uuid='', v_mag=0.0, v_phase=0.0, p=0.0, q=0.0, index=0,  bus_type='PQ'):
-		self.index = index
+	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
 		self.uuid = uuid
-		self.power = complex(p, q)
-		self.voltage = v_mag*np.cos(v_phase) + 1j * v_mag*np.sin(v_phase)
+		self.index = index
+		self.baseVoltage = base_voltage
+		self.base_apparent_power = base_apparent_power
 		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):
-		self.r = r
-		self.x = x
+	def __init__(self, r, x, start_node, end_node, base_voltage=1.0, base_apparent_power=1.0):
+		self.baseVoltage = base_voltage
+		self.base_apparent_power = base_apparent_power
+		self.base_impedance = base_voltage**2 / self.base_apparent_power
 		self.start_node = start_node
 		self.end_node = end_node
+		self.r = r
+		self.x = x
 		self.z = self.r + 1j*self.x
 		self.y = 1/self.z if (self.z != 0) else float("inf")
+		self.r_pu = r/self.base_impedance
+		self.x_pu = x/self.base_impedance
+		self.z_pu = self.r_pu + 1j*self.x_pu
+		self.y_pu = 1/self.z_pu if (self.z_pu != 0) else float("inf")
 
 class System():	
 	def __init__(self):
@@ -36,60 +48,101 @@ class System():
 		self.bX=[]
 		self.P=[]
 		self.Q=[]
-		self.Ymatrix = np.zeros((1, 1),dtype=np.complex)
+		self.Ymatrix = np.zeros([], 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():
-			if value.__class__.__name__=="TopologicalNode":
-				self.P.append(value.pInjection)
-				self.Q.append(value.qInjection)
+	def load_cim_data(self, res, base_apparent_power):
+		"""
+		To fill the vectors node, branch, bR, bX, P and Q
+		"""
+
+		for uuid, element in res.items():
+			if element.__class__.__name__=="TopologicalNode":
+				vmag = 0.0
+				vphase = 0.0
+				pInj = 0.0
+				qinj = 0.0
+				for uuid2, element2 in res.items():
+					if element2.__class__.__name__=="SvVoltage":
+						if element2.getNodeUUID() == uuid:
+							vmag = element2.v
+							vphase = element2.angle
+							break
+				for uuid2, element2 in res.items():
+					if element2.__class__.__name__=="SvPowerFlow":
+						if element2.getNodeUUID() == uuid:
+							pInj += element2.p
+							qinj += element2.q
+				self.P.append(pInj)
+				self.Q.append(qinj)
 				index=len(self.P)-1
-				self.nodes.append(Node(name=value.name, uuid=value.mRID, p=value.pInjection, q=value.qInjection, index=index))
-		for key, value in res.items():
-			if value.__class__.__name__=="ACLineSegment":
-				length=value.length
+				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))
+
+		for uuid, element in res.items():
+			if element.__class__.__name__=="ACLineSegment":
+				length=element.length
 				if length==0.0:
 					length=1.0
-				self.bR.append(value.r*length)
-				self.bX.append(value.x*length)
-				for i in range(len(self.nodes)):
-					if value.startNodeID==self.nodes[i].uuid:
-						startNode=self.nodes[i]
+				bR = element.r*length
+				bX = element.x*length
+				self.bR.append(bR)
+				self.bX.append(bX)
+				for node in self.nodes:
+					if element.startNodeID==node.uuid:
+						startNode = node
 						break
-				for i in range(len(self.nodes)):
-					if value.endNodeID==self.nodes[i].uuid:
-						endNode=self.nodes[i]
+				for node in self.nodes:
+					if element.endNodeID==node.uuid:
+						endNode=node
 						break
-				self.branches.append(Branch(value.r, value.x, startNode, endNode))	
-			elif value.__class__.__name__=="PowerTransformer":
-				self.bR.append(value.primaryConnection.r)
-				self.bX.append(value.primaryConnection.x)
+				
+				base_voltage = element.BaseVoltage.nominalVoltage
+				self.branches.append(Branch(bR, bX, startNode, endNode, base_voltage, base_apparent_power))
+			
+			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 value.startNodeID==self.nodes[i].uuid:
+					if element.startNodeID==self.nodes[i].uuid:
 						startNode=self.nodes[i]
 						break
 				for i in range(len(self.nodes)):
-					if value.endNodeID==self.nodes[i].uuid:
+					if element.endNodeID==self.nodes[i].uuid:
 						endNode=self.nodes[i]
 						break
-				self.branches.append(Branch(value.primaryConnection.r, value.primaryConnection.x, startNode, endNode))
+
+				#base voltage = high voltage side (=primaryConnection)
+				base_voltage = element.primaryConnection.BaseVoltage.nominalVoltage
+				self.branches.append(Branch(bR, bX, startNode, endNode, base_voltage, base_apparent_power))
+			
 			else:
 				continue
 
-		#determine the impedance matrix
-		#self.Ymatrix_calc(Zb)
+		#calculate impedance matrix
 		self.Ymatrix_calc()
 
-	#def Ymatrix_calc(self, Zb):
-	def Ymatrix_calc(self):
+	def _getNodeType(self, node):
 		"""
-		@param Zb: base value of impedance
+		return the type of a node: PQ, PV or SLACK
+
+		@param node: element of class cimpy.Topology.TopologicalNode
 		"""
+		for terminal in node.Terminal:
+			if terminal.ConductingEquipment.__class__.__name__=="ExternalNetworkInjection":
+				return "SLACK"
+			elif terminal.ConductingEquipment.__class__.__name__=="SynchronousMachine":
+				return "PV"
+		return "PQ"
+
+	def Ymatrix_calc(self):
 		nodes_num = len(self.nodes)
-		self.Ymatrix = np.zeros((nodes_num, nodes_num),dtype=np.complex)
+		self.Ymatrix = np.zeros((nodes_num, nodes_num), dtype=np.complex)
 		self.Adjacencies = [[] for _ in range(nodes_num)]
 		for branch in self.branches:
 			fr = branch.start_node.index
@@ -101,17 +154,18 @@ 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()
 	
 	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]))
-		else:
+		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]))

acs/state_estimation/nv_powerflow_cim.py

@@ -60,7 +60,7 @@ def solve(system):
 				h[m] = np.inner(H[m],State)
 				h[m+1] = np.inner(H[m+1],State)
 			elif type is BusType.PQ:
-				z[m] = (np.real(system.nodes[i].power)*np.real(V[i])+ np.imag(system.nodes[i].power)*np.imag(V[i]))/(np.abs(V[i])**2)
+				z[m] = (np.real(system.nodes[i].power)*np.real(V[i]) + np.imag(system.nodes[i].power)*np.imag(V[i]))/(np.abs(V[i])**2)
 				z[m+1] = (np.real(system.nodes[i].power)*np.imag(V[i]) - np.imag(system.nodes[i].power)*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)
@@ -78,16 +78,11 @@ def solve(system):
 		diff = np.amax(np.absolute(Delta_State))
 		
 		V = State[:nodes_num] + 1j * State[nodes_num:]
-		
 		num_iter = num_iter+1
 		
+	# calculate all the other quantities of the grid
 	powerflow_results = results.Results(system)
 	powerflow_results.load_voltages(V)
-	powerflow_results.calculateI()
-	powerflow_results.calculateIinj()
-	powerflow_results.calculateSinj()
-	powerflow_results.calculateI()
-	powerflow_results.calculateS1()
-	powerflow_results.calculateS2()
-	
+	powerflow_results.calculate_all()
+
 	return powerflow_results, num_iter
\ No newline at end of file

acs/state_estimation/nv_state_estimator_cim.py

@@ -662,10 +662,10 @@ 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
+				m2 = m + nodes_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
+				n2 = n + nodes_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
 
+import sys
+sys.path.append("../../../dataprocessing")
+from villas.dataprocessing.readtools import *
+
 class ResultsNode():
 	def __init__(self, topo_node):		
 		self.topology_node = topo_node		
@@ -32,7 +36,7 @@ class Results():
 		"""
 		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]
+			node.voltage = loadflow_results[node.topology_node.uuid].values[0]
 
 	def load_voltages(self, V):
 		"""
@@ -43,15 +47,6 @@ class Results():
 				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
@@ -131,7 +126,7 @@ class Results():
 				
 	def get_voltages(self):
 		"""
-		get complex Power Injection at nodes
+		get node voltages
 		for a test purpose
 		"""
 		voltages = np.zeros(len(self.nodes), dtype=np.complex_)

examples/95bus_tests/test_nv_powerflow_cim.py

 import math
 import sys
 
-
 sys.path.append("../../acs/state_estimation")
 import network
 

examples/CIGREMV_tests/comparison_dpsim_state_estimator.py

@@ -14,7 +14,7 @@ import measurement
 import nv_state_estimator_cim
 import results
 
-logging.basicConfig(filename='CIGRE.log', level=logging.INFO)
+logging.basicConfig(filename='CIGRE.log', level=logging.INFO, filemode='w')
 
 cim_xml_path = r"..\..\..\cim-grid-data\CIGRE_MV\CIGRE_MV_no_tapchanger_With_LoadFlow_Results"
 cim_xml_files=[cim_xml_path + r"\Rootnet_FULL_NE_06J16h_DI.xml", 
@@ -24,8 +24,6 @@ cim_xml_files=[cim_xml_path + r"\Rootnet_FULL_NE_06J16h_DI.xml",
 
 #read cim files and create new network.Systen object
 res=cimpy.cimread(cim_xml_files)
-cimpy.setNodes(res)
-cimpy.setPowerTransformerEnd(res)
 system = network.System()
 system.load_cim_data(res)
 

examples/CIGREMV_tests/test_nv_powerflow_cim.py

+import sys
+import logging
+
+sys.path.append("../../acs/state_estimation")
+import network
+import nv_powerflow_cim
+
+sys.path.append("../../../cimpy")
+import cimpy
+
+logging.basicConfig(filename='CIGRE.log', level=logging.INFO, filemode='w')
+
+cim_xml_path = r"C:\Users\Martin\Desktop\hiwi\git\cim-grid-data\CIGRE_MV\CIGRE_MV_no_tapchanger_With_LoadFlow_Results"
+cim_xml_files=[cim_xml_path + r"\Rootnet_FULL_NE_06J16h_DI.xml", 
+		   cim_xml_path + r"\Rootnet_FULL_NE_06J16h_EQ.xml",
+		   cim_xml_path + r"\Rootnet_FULL_NE_06J16h_SV.xml",
+		   cim_xml_path + r"\Rootnet_FULL_NE_06J16h_TP.xml"]
+
+#read cim files and create new network.Systen object
+res=cimpy.cimread(cim_xml_files)
+system = network.System()
+system.load_cim_data(res, 20)
+
+#print node voltages
+for node in system.nodes:
+    print('{}={}'.format(node.uuid, node.voltage))
+
+print()
+
+#print node powers
+for node in system.nodes:
+    print('{}={}'.format(node.uuid, node.power))
+
+results, num_iter_cim = nv_powerflow_cim.solve(system)
+
+print()
+print("voltages:")
+#print node voltages
+for node in results.nodes:
+    print('{}={}'.format(node.topology_node.uuid, node.voltage))
+
+#results.print_voltages_polar()
\ No newline at end of file