improve state estimation and fixed some errors

acs/state_estimation/nv_powerflow_cim.py

@@ -98,7 +98,15 @@ class PowerflowResults():
 			for node in self.nodes:
 				if branch_index == node.topology_node.index:
 					branch.power2 = -node.voltage*(np.conj(branch.current))
-					
+	
+	def get_node(self, index):
+		"""
+		return the PowerflowNode with PowerflowNode.topology_node.index == index
+		"""
+		for node in self.nodes:
+			if index == node.topology_node.index:
+				return node
+				
 	def get_voltages(self):
 		"""
 		get complex Power Injection at nodes

acs/state_estimation/nv_state_estimator_cim.py

 import sys
 import numpy
+from enum import Enum
 
 sys.path.append("../../../dataprocessing")
 from villas.dataprocessing.readtools import read_timeseries_dpsim
+import nv_powerflow_cim
+
+class ElemType(Enum):
+	Node = 1		#Node Voltage
+	Branch = 2		#Complex Power Injection at node
+	
+class MeasType(Enum):
+	V = 1		#Node Voltage
+	Sinj= 2		#Complex Power Injection at node
+	S1 = 4		#Complex Power flow at branch, measured at initial node
+	S2 = 5		#Complex Power flow at branch, measured at final node
+	I = 6		#Branch Current
+	Vpmu = 7	#Node Voltage
+	Ipmu = 9	#Branch Current
 
 class Measurement():
-	def __init__(self, topo_node, meas_type, meas_value, std_dev):
+	def __init__(self, element, element_type, 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: 
+		@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.
 		"""
 		
-		if meas_type not in ("v", "Vpmu", "p", "q", "i", "Ipmu"):
-			raise Exception("Invalid measurement type (%s)" % meas_type)
+		if not isinstance(element_type, ElemType):
+			raise Exception("elem_type must be an object of class ElemType")
+		
+		if not isinstance(meas_type, MeasType):
+			raise Exception("meas_type must be an object of class MeasType")
 			
-		self.topology_node = topo_node
+		self.element = element
 		self.meas_type = meas_type
+		self.element_type = element_type
 		self.meas_value = meas_value
 		self.std_dev = self.meas_value*(std_dev/300)
 
@@ -26,8 +46,8 @@ 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 create_measurement(self, element, element_type, meas_type, meas_value, std_dev):
+		self.measurements_set.append(Measurement(element, element_type, meas_type, meas_value, std_dev))
 	
 def DsseCall(system, zdata, measurements):
 	""" It identifies the type of measurements present in the measurement set and 
@@ -37,7 +57,6 @@ def DsseCall(system, zdata, measurements):
 
 	# 1. Select type of Estimator.
 	# If at least a PMU is present we launch the combined estimator, otherwise a simple traditional etimator 
-
 	Vmag_meas = 0
 	Vpmu_meas = 0
 	for elem in measurements.measurements_set:
@@ -48,12 +67,9 @@ def DsseCall(system, zdata, measurements):
 	
 	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)
 	elif est_code == 2:
@@ -70,48 +86,16 @@ def DsseCall(system, zdata, measurements):
 	# S2est:  Complex Power flow at branch, measured at final node  
 	# Sinjest: Complex Power Injection at node
 
-	nodes_num = len(system.nodes)
-	branches_num = len(system.branches)
-	
 	""" From here on, all the other quantities of the grid are calculated """
-	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])*system.Ymatrix[fr][to]
-	Ir = numpy.real(Irx)
-	Ix = numpy.imag(Irx)
-	
-	#Iest = Real_to_all(Ir,Ix)
-	Iinj_r = numpy.zeros(nodes_num)
-	Iinj_x = numpy.zeros(nodes_num)
-	for k in range(nodes_num):
-		to=[]
-		fr=[]
-		for m in range(branches_num):
-			if k==system.branches[m].start_node.index:
-				fr.append(m)
-			if k==system.branches[m].end_node.index:
-				to.append(m)
-		
-		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)
-	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))
-	
-	S1est=numpy.array([])
-	S2est=numpy.array([])
-	for i in range(branches_num):
-		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))
-		
-	return Vest, Iest, Iinjest, S1est, S2est, Sinjest
+	results = nv_powerflow_cim.PowerflowResults(system)
+	results.load_voltages(Vest)
+	results.calculateI()
+	results.calculateIinj()
+	results.calculateSinj()
+	results.calculateI()
+	results.calculateS()
+	
+	return results
 
 def DsseTrad(system, zdata):
 	""" It performs state estimation using rectangular node voltage state variables 
@@ -240,7 +224,7 @@ def DsseTrad(system, zdata):
 		
 		""" Voltage Magnitude Measurements """
 		# at every iteration we update h(x) vector where V measure are available
-		h1 = V[busvi-1].mag[busvi-1]
+		h1 = numpy.absolute(V[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):

examples/95bus_tests/test_nv_powerflow_cim.py

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

examples/95bus_tests/test_nv_state_estimation_cim_new.py

+import numpy as np
+import math
+import sys
+import copy 
+
+sys.path.append("../../acs/state_estimation")
+import py_95bus_network_data
+import py_95bus_meas_data
+import nv_state_estimator
+import nv_powerflow
+
+import network
+import nv_powerflow_cim
+import nv_state_estimator_cim
+
+class PerUnit:
+    def __init__(self, S, V):
+        self.S = S
+        self.V = V
+        self.I = S/V
+        self.Z = S/(V**2)
+        
+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)
+		
+		
+""" Insert here per unit values of the grid for power and voltage """
+S = 100*(10**6)
+V = (11*(10**3))/math.sqrt(3)
+slackV = 1.02
+
+Base = PerUnit(S,V)
+branch, node = py_95bus_network_data.Network_95_nodes(Base, slackV)
+Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue, num_iter = nv_powerflow.solve(branch, node)
+
+system = network.load_python_data(node, branch, node.type)
+res, num_iter = nv_powerflow_cim.solve(system)
+
+""" Write here the indexes of the nodes/branches where the measurements are"""
+V_idx = np.array([1,11,55])
+I_idx = np.array([13,36])
+Sinj_idx = np.linspace(2,node.num,node.num-1)
+S1_idx = np.array([1,11,28,55,59])
+S2_idx = np.array([10,54])
+Ipmu_idx = np.array([13,36])
+Vpmu_idx = np.array([1,11,55])
+
+""" Write here the percent uncertainties of the measurements""" 
+V_unc = 1
+I_unc = 2
+Sinj_unc = 5
+S_unc = 2
+Pmu_mag_unc = 0.7
+Pmu_phase_unc = 0.7
+
+""" create measurements set (for nv_state_estimator_cim.py)"""
+measurements_set = nv_state_estimator_cim.Measurents_set()
+for i in [0,10,54]:
+	PowerflowNode = res.get_node(i)
+	measurements_set.create_measurement(PowerflowNode.topology_node, nv_state_estimator_cim.ElemType.Node, nv_state_estimator_cim.MeasType.V, PowerflowNode.voltage, V_unc)
+	measurements_set.create_measurement(PowerflowNode.topology_node, nv_state_estimator_cim.ElemType.Node, nv_state_estimator_cim.MeasType.Vpmu, PowerflowNode.voltage, Pmu_mag_unc)
+for i in [12,35]:
+	measurements_set.create_measurement(res.branches[i].topology_branch, nv_state_estimator_cim.ElemType.Branch,  nv_state_estimator_cim.MeasType.I, res.branches[i].current, I_unc)
+	measurements_set.create_measurement(res.branches[i].topology_branch, nv_state_estimator_cim.ElemType.Branch,  nv_state_estimator_cim.MeasType.Ipmu, res.branches[i].current, Pmu_mag_unc)
+for i in range(1,len(res.nodes)):
+	PowerflowNode = res.get_node(i)
+	measurements_set.create_measurement(PowerflowNode.topology_node, nv_state_estimator_cim.ElemType.Node, nv_state_estimator_cim.MeasType.Sinj, PowerflowNode.power, Sinj_unc)
+for i in [0,10,27,54,58]:
+	measurements_set.create_measurement(res.branches[i].topology_branch, nv_state_estimator_cim.ElemType.Branch,  nv_state_estimator_cim.MeasType.S1, res.branches[i].power, I_unc)
+for i in [9,53]:
+	measurements_set.create_measurement(res.branches[i].topology_branch, nv_state_estimator_cim.ElemType.Branch,  nv_state_estimator_cim.MeasType.S2, res.branches[i].power2, I_unc)
+
+print("########################   TEST   ####################")
+Meas_mag = np.array([])
+Meas_phase = np.array([])
+for meas in measurements_set.measurements_set:
+	if meas.meas_type == nv_state_estimator_cim.MeasType.Vpmu:
+		Meas_mag = np.append(Meas_mag, np.absolute(meas.meas_value))
+		Meas_phase = np.append(Meas_phase, np.angle(meas.meas_value))
+		#Meas_mag = np.append(Meas_mag, meas.meas_value.real)
+		#Meas_phase = np.append(Meas_phase, meas.meas_value.imag)
+print(Meas_mag)
+print(Meas_phase)
+
+	
+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)
+zdatameas =  Zdata_init(meas)
+
+zdata, zdatameas = py_95bus_meas_data.Zdatatrue_creation(zdata, zdatameas, meas, branch, Vtrue, Itrue, Iinjtrue, S1true, S2true, Sinjtrue)
\ No newline at end of file