restructure files and refactoring of state estimator and dpsim comparison

.gitignore

@@ -4,4 +4,7 @@
 
 # python folders
 *.egg-info/
-__pycache__
\ No newline at end of file
+__pycache__
+
+# log files
+*.log
\ No newline at end of file

acs/state_estimation/measurement_generator.py

+import numpy as np
+
+class Measurements:
+	def __init__(self, index, unc):
+		self.index = index.astype(int)
+		self.unc = unc
+		self.num = len(index)
+		
+class Measurement_set:
+	def __init__(self, V, I, Sinj, S1, S2, Vpmu_mag, Vpmu_phase, Ipmu_mag, Ipmu_phase):
+		self.V = V
+		self.Sinj = Sinj
+		self.S1 = S1
+		self.S2 = S2
+		self.I = I
+		self.Vpmu_mag = Vpmu_mag
+		self.Vpmu_phase = Vpmu_phase 
+		self.Ipmu_mag = Ipmu_mag
+		self.Ipmu_phase = Ipmu_phase
+		
+class Zdata_init:
+	def __init__(self, meas):
+		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 = np.zeros(nmeas)
+		self.mval = np.zeros(nmeas)
+		self.mbranch = np.zeros(nmeas)
+		self.mfrom = np.zeros(nmeas)
+		self.mto = np.zeros(nmeas)
+		self.mstddev = np.zeros(nmeas)
+
+def Zdatatrue_creation(zdata, zdatameas, meas, system, V, I, Iinj, S1, S2, Sinj):
+	""" It gives the true values at the measurement points."""
+	import numpy as np
+	
+	t1 = np.ones(meas.V.num)
+	t2 = 2*np.ones(meas.Sinj.num)
+	t3 = 3*np.ones(meas.Sinj.num)
+	t4 = 4*np.ones(meas.S1.num + meas.S2.num)
+	t5 = 5*np.ones(meas.S1.num + meas.S2.num)
+	t6 = 6*np.ones(meas.I.num)
+	t7 = 7*np.ones(meas.Vpmu_mag.num)
+	t8 = 8*np.ones(meas.Vpmu_phase.num)
+	t9 = 9*np.ones(meas.Ipmu_mag.num)
+	t10 = 10*np.ones(meas.Ipmu_phase.num)
+	
+	zdata.mtype = np.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
+	zdatameas.mtype = np.concatenate((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10),axis=0)
+	
+	z1 = np.abs(V[meas.V.index-1])
+	z2 = Sinj.real[meas.Sinj.index-1]
+	z3 = Sinj.imag[meas.Sinj.index-1]
+	z4_1 = S1.real[meas.S1.index-1]
+	z4_2 = S2.real[meas.S2.index-1]
+	z5_1 = S1.imag[meas.S1.index-1]
+	z5_2 = S2.imag[meas.S2.index-1]
+	z6 = np.abs(I[meas.I.index-1])
+	z7 = np.abs(V[meas.Vpmu_mag.index-1])
+	z8 = np.angle(V[meas.Vpmu_phase.index-1])
+	z9 = np.abs(I[meas.Ipmu_mag.index-1])
+	z10 = np.angle(I[meas.Ipmu_phase.index-1])
+	
+	zdata.mval = np.concatenate((z1,z2,z3,z4_1,z4_2,z5_1,z5_2,z6,z7,z8,z9,z10),axis=0)
+	
+	"""
+	b1 = np.zeros(meas.V.num)
+	b2 = np.zeros(meas.Sinj.num)
+	b3 = np.zeros(meas.Sinj.num)
+	b4_1 = branch.code[meas.S1.index-1]
+	b4_2 = -branch.code[meas.S2.index-1]
+	b5_1 = branch.code[meas.S1.index-1]
+	b5_2 = -branch.code[meas.S2.index-1]
+	b6 = branch.code[meas.I.index-1]
+	b7 = np.zeros(meas.Vpmu_mag.num)
+	b8 = np.zeros(meas.Vpmu_phase.num)
+	b9 = branch.code[meas.Ipmu_mag.index-1]
+	b10 = branch.code[meas.Ipmu_phase.index-1]
+	
+	zdata.mbranch = np.concatenate((b1,b2,b3,b4_1,b4_2,b5_1,b5_2,b6,b7,b8,b9,b10),axis=0)
+	zdata.mbranch = zdata.mbranch.astype(int)
+	zdatameas.mbranch = np.concatenate((b1,b2,b3,b4_1,b4_2,b5_1,b5_2,b6,b7,b8,b9,b10),axis=0)
+	zdatameas.mbranch = zdata.mbranch.astype(int)
+	"""
+	
+	fr1 = meas.V.index
+	fr2 = meas.Sinj.index
+	fr3 = meas.Sinj.index
+	fr4_1 = np.array([])
+	fr4_2 = np.array([])
+	fr5_1 = np.array([])
+	fr5_2 = np.array([])
+	for index in range(len(meas.S1.index)):
+		fr4_1 = np.append(fr4_1, system.branches[meas.S1.index[index]-1].start_node.index)
+		fr5_1 = np.append(fr5_1, system.branches[meas.S1.index[index]-1].start_node.index)
+	for index in range(len(meas.S2.index)):
+		fr4_2 = np.append(fr4_2, system.branches[meas.S2.index[index]-1].end_node.index)
+		fr5_2 = np.append(fr5_2, system.branches[meas.S2.index[index]-1].end_node.index)
+	fr6 = np.array([])
+	for index in range(len(meas.I.index)):
+		fr6 = np.append(fr6, system.branches[meas.I.index[index]-1].start_node.index)
+	fr7 = meas.Vpmu_mag.index
+	fr8 = meas.Vpmu_phase.index
+	fr9 = np.array([])
+	fr10 = np.array([])
+	for index in range(len(meas.Ipmu_mag.index)):
+		fr9 = np.append(fr9, system.branches[meas.Ipmu_mag.index[index]-1].start_node.index)
+	for index in range(len(meas.Ipmu_phase.index)):
+		fr10 = np.append(fr10, system.branches[meas.Ipmu_phase.index[index]-1].start_node.index)	
+	
+	zdata.mfrom = np.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
+	zdata.mfrom = zdata.mfrom.astype(int)
+	zdatameas.mfrom = np.concatenate((fr1,fr2,fr3,fr4_1,fr4_2,fr5_1,fr5_2,fr6,fr7,fr8,fr9,fr10),axis=0)
+	zdatameas.mfrom = zdata.mfrom.astype(int)
+	
+	to1 = np.zeros(meas.V.num)
+	to2 = np.zeros(meas.Sinj.num)
+	to3 = np.zeros(meas.Sinj.num)
+	to4_1 = np.array([])
+	to4_2 = np.array([])
+	to5_1 = np.array([])
+	to5_2 = np.array([])
+	for index in range(len(meas.S1.index)):
+		to4_1 = np.append(to4_1, system.branches[meas.S1.index[index]-1].end_node.index)
+		to5_1 = np.append(to5_1, system.branches[meas.S1.index[index]-1].end_node.index)
+	for index in range(len(meas.S2.index)):
+		to4_2 = np.append(to4_2, system.branches[meas.S2.index[index]-1].start_node.index)
+		to5_2 = np.append(to5_2, system.branches[meas.S2.index[index]-1].start_node.index)
+	to6 = np.array([])
+	for index in range(len(meas.I.index)):
+		to6 = np.append(to6, system.branches[meas.I.index[index]-1].end_node.index)
+	to7 = np.zeros(meas.Vpmu_mag.num)
+	to8 = np.zeros(meas.Vpmu_phase.num)
+	to9 = np.array([])
+	to10 = np.array([])
+	for index in range(len(meas.Ipmu_mag.index)):
+		to9 = np.append(to9, system.branches[meas.Ipmu_mag.index[index]-1].end_node.index)
+	for index in range(len(meas.Ipmu_phase.index)):
+		to10 = np.append(to10, system.branches[meas.Ipmu_phase.index[index]-1].end_node.index)
+	
+	zdata.mto = np.concatenate((to1,to2,to3,to4_1,to4_2,to5_1,to5_2,to6,to7,to8,to9,to10),axis=0)
+	zdata.mto = zdata.mto.astype(int)
+	zdatameas.mto = np.concatenate((to1,to2,to3,to4_1,to4_2,to5_1,to5_2,to6,to7,to8,to9,to10),axis=0)
+	zdatameas.mto = zdata.mto.astype(int)
+	
+	d1 = z1*(meas.V.unc/300)
+	d2 = np.absolute(z2*(meas.Sinj.unc/300))
+	d3 = np.absolute(z3*(meas.Sinj.unc/300))
+	d4_1 = np.absolute(z4_1*(meas.S1.unc/300))
+	d4_2 = np.absolute(z4_2*(meas.S2.unc/300))
+	d5_1 = np.absolute(z5_1*(meas.S1.unc/300))
+	d5_2 = np.absolute(z5_2*(meas.S2.unc/300))
+	d6 = z6*(meas.I.unc/300)
+	d7 = z7*(meas.Vpmu_mag.unc/300)
+	d8 = (meas.Vpmu_phase.unc/300)*np.ones(meas.Vpmu_phase.num)
+	d9 = z9*(meas.Ipmu_mag.unc/300)
+	d10 = (meas.Ipmu_phase.unc/300)*np.ones(meas.Ipmu_phase.num)
+	
+	zdata.mstddev = np.concatenate((d1,d2,d3,d4_1,d4_2,d5_1,d5_2,d6,d7,d8,d9,d10),axis=0)
+	zdatameas.mstddev = np.concatenate((d1,d2,d3,d4_1,d4_2,d5_1,d5_2,d6,d7,d8,d9,d10),axis=0)
+#	idx = np.where(zdata.mstddev<10**(-6))
+#	zdata.mstddev[idx] = 10**(-6)
+	
+	return zdata, zdatameas
+
+def Zdatameas_creation(zdata, zdatameas):
+	""" It gives the measured values (affected by uncertainty) at the measurement points."""
+	import numpy as np
+	
+	zmeas = zdata.mval
+	zdev = zdata.mstddev
+	err_pu = np.random.normal(0,1,len(zmeas))
+	zdatameas.mval = zmeas + np.multiply(zdev,err_pu)
+	
+	return zdatameas
\ No newline at end of file

examples/CIGREMV_tests/comparison_dpsim_state_estimator.py

+import numpy as np
+import math
+import sys
+import copy
+import logging
+import matplotlib.pyplot as plt
+
+sys.path.append("../../cimpy")
+import cimpy
+
+sys.path.append("./acs/state_estimation")
+import network
+import nv_powerflow_cim
+from measurement_generator import *
+from dpsim_reader import * # TODO replace with functions from villas.dataprocessing
+import nv_state_estimator_cim
+
+logging.basicConfig(filename='CIGRE.log', level=logging.INFO)
+
+cim_xml_path = r"D:\git\data\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"]
+
+loadflow_results_path = r"D:\git\data\reference-results\DPsim\StaticPhasor"
+loadflow_results_file = loadflow_results_path + r"\CIGRE-MV-NoTap.csv" 
+
+res=cimpy.cimread(cim_xml_files)
+cimpy.setNodes(res)
+cimpy.setPowerTransformerEnd(res)
+system = network.System()
+system.load_cim_data(res)
+Ymatr, Adj = network.Ymatrix_calc(system)
+
+loadflow_results = read_data(loadflow_results_file)
+
+#order readed data according to system.nodes
+Vtrue=np.zeros(len(loadflow_results), dtype=np.complex_)
+for elem in range(len(system.nodes)): 
+	Vtrue[elem] = loadflow_results[system.nodes[elem].uuid]
+
+""" Write here the indexes of the nodes/branches where the measurements are"""
+V_idx = np.array([])
+I_idx = np.array([])
+Sinj_idx = np.array([])
+S1_idx = np.array([])
+S2_idx = np.array([])
+Ipmu_idx = np.array([])
+Vpmu_idx= np.linspace(1, len(Vtrue), len(Vtrue))
+
+# --- State Estimation with Ideal Measurements ---
+""" Write here the percent uncertainties of the measurements""" 
+V_unc = 0
+I_unc = 0
+Sinj_unc = 0
+S_unc = 0
+Pmu_mag_unc = 0
+Pmu_phase_unc = 0
+
+# Create measurements data structures
+V = Measurements(V_idx,V_unc)
+I = Measurements(I_idx,I_unc)
+Sinj = Measurements(Sinj_idx,Sinj_unc)
+S1 = Measurements(S1_idx,S_unc)
+S2 = Measurements(S2_idx,S_unc)
+Ipmu_mag = Measurements(Ipmu_idx,Pmu_mag_unc)
+Ipmu_phase = Measurements(Ipmu_idx,Pmu_phase_unc)
+Vpmu_mag = Measurements(Vpmu_idx,Pmu_mag_unc)
+Vpmu_phase = Measurements(Vpmu_idx,Pmu_phase_unc)
+
+meas_ideal = Measurement_set(V, I, Sinj, S1, S2, Vpmu_mag, Vpmu_phase, Ipmu_mag, Ipmu_phase)
+zdata_ideal = Zdata_init(meas_ideal)
+zdata_ideal.mtype = np.ones(meas_ideal.V.num)
+zdatameas_ideal =  Zdata_init(meas_ideal)
+
+zdata_ideal, zdatameas_ideal = Zdatatrue_creation(zdata_ideal, zdatameas_ideal, meas_ideal, system, Vtrue, np.array([]), np.array([]), np.array([]), np.array([]), np.array([]))
+zdatameas_ideal = Zdatameas_creation(zdata_ideal, zdatameas_ideal)
+
+# Perform state estimation
+Vest_ideal, Iest, Iinjest, S1est, S2est, Sinjest = nv_state_estimator_cim.DsseCall(system, zdatameas_ideal, Ymatr, Adj)
+
+# Show numerical comparison
+print(Vest_ideal.complex - Vtrue)
+
+# --- State Estimation with Non-Ideal Measurements ---
+""" Write here the percent uncertainties of the measurements""" 
+Pmu_mag_unc = 1
+
+# Create measurements data structures
+Vpmu_mag = Measurements(Vpmu_idx,Pmu_mag_unc)
+
+meas_real = Measurement_set(V, I, Sinj, S1, S2, Vpmu_mag, Vpmu_phase, Ipmu_mag, Ipmu_phase)
+zdata_real = Zdata_init(meas_real)
+zdata_real.mtype = np.ones(meas_real.V.num)
+zdatameas_real =  Zdata_init(meas_real)
+
+zdata_real, zdatameas_real = Zdatatrue_creation(zdata_real, zdatameas_real, meas_real, system, Vtrue, np.array([]), np.array([]), np.array([]), np.array([]), np.array([]))
+zdatameas_real = Zdatameas_creation(zdata_real, zdatameas_real)
+
+# Perform state estimation
+Vest_real, Iest, Iinjest, S1est, S2est, Sinjest = nv_state_estimator_cim.DsseCall(system, zdatameas_real, Ymatr, Adj)
+
+# Show numerical comparison
+print(Vest_real.complex - Vtrue)
+
+# Plot comparison
+line_width = 3
+fontsize = 22
+
+plt.figure()
+ax = plt.subplot()
+
+nodes_uuid = [system.nodes[elem].uuid for elem in range(len(system.nodes))]
+
+# Reorder and rescale results
+idx_filter = np.argsort([int(uuid[1:]) for uuid in nodes_uuid])[1:]
+nodes_uuid_filtered = [nodes_uuid[idx] for idx in idx_filter]
+Vtrue_filtered = [abs(Vtrue[idx]/1000) for idx in idx_filter]
+Vest_ideal_filtered = [abs(Vest_ideal.complex[idx]/1000) for idx in idx_filter]
+Vest_real_filtered = [abs(Vest_real.complex[idx]/1000) for idx in idx_filter]
+
+
+plt.plot(Vest_ideal_filtered, linewidth=line_width, linestyle='-', label="State estimator (ideal measurements)")
+plt.plot(Vtrue_filtered, linewidth=line_width, linestyle=':', label="DPsim load flow results")
+plt.plot(Vest_real_filtered, linewidth=line_width, linestyle='-', label="State estimator (real measurements)")
+
+plt.xticks(range(len(system.nodes)), fontsize=fontsize)
+plt.yticks(fontsize=fontsize)
+ax.set_xticklabels(nodes_uuid_filtered)
+plt.ylabel("Node Voltage [kV]", fontsize=fontsize)
+plt.xlim([0,len(system.nodes)-2])
+plt.legend(fontsize=fontsize-2)
+plt.show()
\ No newline at end of file

examples/CIGREMV_tests/dpsim_reader.py

+def read_data(file_name):
+	file = []
+	with open(file_name, "r") as filestream:
+		for line in filestream:
+			file.append([x.strip() for x in line.split(',')])
+	
+	#remove "time"
+	file[0].pop(0) 
+	file[1].pop(0)
+	
+	results={}
+	for elem in range(len(file[0])):
+		node, type = file[0][elem].split('.')
+		if not node in results:
+			results[node] = 0.0 + 0.0j
+		if type=="real":
+			results[node] += float(file[1][elem])
+		elif type=="imag":
+			results[node] += complex(0,float(file[1][elem]))
+	
+	return results
\ No newline at end of file

examples/cim_CIGRE_state_estimation_test.py

-import numpy as np
-import math
-import sys
-import copy
-import logging
-
-sys.path.append("../../cimpy")
-import cimpy
-
-sys.path.append("../acs/state_estimation")
-import network
-import nv_powerflow_cim
-import cim_py_95bus_meas_data
-import nv_state_estimator_cim
-
-logging.basicConfig(filename='CIGRE.log', level=logging.INFO)
-
-class Measurements:
-	def __init__(self, index, unc):
-		self.index = index.astype(int)
-		self.unc = unc
-		self.num = len(index)
-		
-class Measurement_set:
-	def __init__(self, V, I, Sinj, S1, S2, Vpmu_mag, Vpmu_phase, Ipmu_mag, Ipmu_phase):
-		self.V = V
-		self.Sinj = Sinj
-		self.S1 = S1
-		self.S2 = S2
-		self.I = I
-		self.Vpmu_mag = Vpmu_mag
-		self.Vpmu_phase = Vpmu_phase 
-		self.Ipmu_mag = Ipmu_mag
-		self.Ipmu_phase = Ipmu_phase
-		
-class Zdata_init:
-	def __init__(self, meas):
-		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 = np.zeros(nmeas)
-		self.mval = np.zeros(nmeas)
-		self.mbranch = np.zeros(nmeas)
-		self.mfrom = np.zeros(nmeas)
-		self.mto = np.zeros(nmeas)
-		self.mstddev = np.zeros(nmeas)
-
-xml_files=[r"..\..\cim-grid-data\CIGRE_MV\CIGRE_MV_no_tapchanger_With_LoadFlow_Results\Rootnet_FULL_NE_06J16h_DI.xml", 
-		   r"..\..\cim-grid-data\CIGRE_MV\CIGRE_MV_no_tapchanger_With_LoadFlow_Results\Rootnet_FULL_NE_06J16h_EQ.xml",
-		   r"..\..\cim-grid-data\CIGRE_MV\CIGRE_MV_no_tapchanger_With_LoadFlow_Results\Rootnet_FULL_NE_06J16h_SV.xml",
-		   r"..\..\cim-grid-data\CIGRE_MV\CIGRE_MV_no_tapchanger_With_LoadFlow_Results\Rootnet_FULL_NE_06J16h_TP.xml"]
-
-res=cimpy.cimread(xml_files)
-cimpy.setNodes(res)
-cimpy.setPowerTransformerEnd(res)
-system = network.System()
-system.load_cim_data(res)
-
-def read_data(file_name):
-	file = []
-	with open(file_name, "r") as filestream:
-		for line in filestream:
-			file.append([x.strip() for x in line.split(',')])
-	
-	#remove "time"
-	file[0].pop(0) 
-	file[1].pop(0)
-	
-	results={}
-	for elem in range(len(file[0])):
-		node, type = file[0][elem].split('.')
-		if not node in results:
-			results[node] = 0.0 + 0.0j
-		if type=="real":
-			results[node] += float(file[1][elem])
-		elif type=="imag":
-			results[node] += complex(0,float(file[1][elem]))
-	
-	return results
-	
-file = r"..\..\reference-results\DPsim\StaticPhasor\CIGRE-MV-NoTap.csv" 
-measurements = read_data(file)
-
-#order readed data according to system.nodes
-Vtrue=np.zeros(len(measurements), dtype=np.complex_)
-for elem in range(len(system.nodes)): 
-	Vtrue[elem] = measurements[system.nodes[elem].uuid]
-
-""" Write here the indexes of the nodes/branches where the measurements are"""
-V_idx = np.linspace(1, len(measurements), len(measurements))
-I_idx = np.array([])
-Sinj_idx = np.array([])
-S1_idx = np.array([])
-S2_idx = np.array([])
-Ipmu_idx = np.array([])
-#Vpmu_idx = np.array([])
-Vpmu_idx=copy.deepcopy(V_idx)
-
-""" Write here the percent uncertainties of the measurements""" 
-V_unc = 0
-I_unc = 0
-Sinj_unc = 0
-S_unc = 0
-Pmu_mag_unc = 0
-Pmu_phase_unc = 0
-
-V = Measurements(V_idx,V_unc)
-I = Measurements(I_idx,I_unc)
-Sinj = Measurements(Sinj_idx,Sinj_unc)
-S1 = Measurements(S1_idx,S_unc)
-S2 = Measurements(S2_idx,S_unc)
-Ipmu_mag = Measurements(Ipmu_idx,Pmu_mag_unc)
-Ipmu_phase = Measurements(Ipmu_idx,Pmu_phase_unc)
-Vpmu_mag = Measurements(Vpmu_idx,Pmu_mag_unc)
-Vpmu_phase = Measurements(Vpmu_idx,Pmu_phase_unc)
-meas = Measurement_set(V, I, Sinj, S1, S2, Vpmu_mag, Vpmu_phase, Ipmu_mag, Ipmu_phase)
-zdata = Zdata_init(meas)
-zdata.mtype = np.ones(meas.V.num)
-zdatameas =  Zdata_init(meas)
-
-zdata, zdatameas = cim_py_95bus_meas_data.Zdatatrue_creation(zdata, zdatameas, meas, system, Vtrue, np.array([]), np.array([]), np.array([]), np.array([]), np.array([]))
-
-MC_trials = 1
-iter_counter = np.zeros(MC_trials)
-Ymatr_cim, Adj_cim = network.Ymatrix_calc(system)
-
-for mciter in range(0,MC_trials):	  
-	zdatameas = cim_py_95bus_meas_data.Zdatameas_creation(zdata, zdatameas)
-	Vest_cim, Iest_cim, Iinjest_cim, S1est_cim, S2est_cim, Sinjest_cim = nv_state_estimator_cim.DsseCall(system, zdatameas, Ymatr_cim, Adj_cim)
-
-print(Vest_cim.complex - Vtrue)
\ No newline at end of file

setup.py

 #!/usr/bin/env python
 
-from distutils.core import setup
+from setuptools import setup
 
 setup(name='acs-state-estimation',
       version='0.1',