Merge branch 'mmo' into 'master'

Mmo

See merge request acs/core/simulation/state-estimation!6

acs/state_estimation/nv_powerflow_cim.py

 import numpy as np
 import math
+from network import Ymatrix_calc
 from network import BusType
 
-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
+class Results():
+	def __init__(self):
+		self.V = []
+		self.I = []
+		self.Iinj = [] 
+		self.S1 = []
+		self.S2 = []
+		self.SInj = []
+		self.num_iter = []
+		
 
 def solve(system):
 	"""It performs Power Flow by using rectangular node voltage state variables."""
@@ -132,4 +128,56 @@ def solve(system):
 		S1=np.append(S1, V[system.branches[i].start_node.index]*(I_conj[i]))
 		S2=np.append(S2, -V[system.branches[i].end_node.index]*(I_conj[i]))
 
-	return V, I, Iinj, S1, S2, Sinj, num_iter
\ No newline at end of file
+	return V, I, Iinj, S1, S2, Sinj, num_iter
+	
+def calculateI(network, V):
+	"""
+	To calculate the branch currents
+	"""
+	I = np.zeros((len(system.branches)), dtype=np.complex)
+	for idx in range(branches_num):
+		fr = system.branches[idx].start_node.index
+		to = system.branches[idx].end_node.index
+		I[idx] = - (V[fr] - V[to])*Ymatrix[fr][to]
+
+def calculateInj(network, I):
+	"""
+	To calculate current injections at a node
+	"""
+	Iinj = np.zeros((len(system.nodes)), dtype=np.complex)
+	for k in range(0, (len(system.nodes))):
+		to=[]
+		fr=[]
+		for m in range(len(system.branches)):
+			if k==system.branches[m].start_node.index:
+				fr.append(m)
+			if k==system.branches[m].end_node.index:
+				to.append(m)
+		Iinj = np.sum(I[to]) - np.sum(I[fr])
+	
+	return Iinj
+
+def calculateS1(system, V, I):
+	"""
+	To calculate powerflow on branches
+	"""
+	S1 = np.zeros((len(system.branches)), dtype=np.complex)
+	for i in range(0, len(system.branches)):
+		S1 = np.append(S1, V[system.branches[i].start_node.index]*(np.conj(I[i])))
+	return S1
+	
+def calculateS2():
+	"""
+	To calculate powerflow on branches
+	"""
+	S2 = np.zeros((len(system.branches)), dtype=np.complex)
+	for i in range(0, len(system.branches)):
+		S2 = np.append(S2, -V[system.branches[i].end_node.index]*(I_conj[i]))
+	return S2
+	
+def calculateSinj(V, Iinj):
+	"""
+	To calculate power injection at a node
+	"""
+	Sinj = np.multiply(V, np.conj(Iinj))
+	return Sinj
\ No newline at end of file

acs/state_estimation/nv_powerflow_cim_v2.py

+import sys
+import math
+import numpy as np
+from network import Ymatrix_calc
+from network import BusType
+
+sys.path.append("../../../dataprocessing")
+from dpsim_reader import read_data
+
+class PF_Results():
+	def __init__(self):
+		self.V = []
+		self.I = []
+		self.Iinj = [] 
+		self.S1 = []
+		self.S2 = []
+		self.SInj = []
+		self.num_iter = 0
+		
+	def read_data(self, file_name, system):
+		"""
+		To read the voltages from a input file
+		"""
+		loadflow_results = read_data(loadflow_results_file)
+		#order readed data according to system.nodes
+		self.V = np.zeros(len(loadflow_results), dtype=np.complex_)
+		for elem in range(len(system.nodes)): 
+			self.V[elem] = loadflow_results[system.nodes[elem].uuid]
+
+def solve(system):
+	"""It performs Power Flow by using rectangular node voltage state variables."""
+	
+	Ymatrix, Adj = Ymatrix_calc(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)))
+	
+	for i in range(0,nodes_num):
+		m = 2*i
+		i2 = i + nodes_num
+		type = system.nodes[i].type 
+		if  type is BusType.SLACK:
+			z[m] = np.real(system.nodes[i].voltage)
+			z[m+1] = np.imag(system.nodes[i].voltage)
+			H[m][i] = 1
+			H[m+1][i2] = 1
+		elif type is BusType.PQ:
+			H[m][i] = - np.real(Ymatrix[i][i])
+			H[m][i2] = np.imag(Ymatrix[i][i])
+			H[m+1][i] = - np.imag(Ymatrix[i][i])
+			H[m+1][i2] = - np.real(Ymatrix[i][i])
+			idx1 = np.subtract(Adj[i],1)
+			idx2 = idx1 + nodes_num
+			H[m][idx1] = - np.real(Ymatrix[i][idx1])
+			H[m][idx2] = np.imag(Ymatrix[i][idx1])
+			H[m+1][idx1] = - np.imag(Ymatrix[i][idx1])
+			H[m+1][idx2] = - np.real(Ymatrix[i][idx1])
+		elif type is BusType.PV:
+			z[m+1] = np.real(system.nodes[i].power)
+			H[m][i] = - np.real(Ymatrix[i][i])
+			H[m][i2] = np.imag(Ymatrix[i][i])
+			idx1 = np.subtract(Adj[i],1)
+			idx2 = idx1 + nodes_num
+			H[m][idx1] = - np.real(Ymatrix[i][idx1])
+			H[m][idx2] = np.imag(Ymatrix[i][idx1])
+	
+	epsilon = 10**(-10)
+	diff = 5
+	V = np.ones(nodes_num) + 1j* np.zeros(nodes_num)
+	num_iter = 0
+	
+	State = np.ones(2*branches_num)
+	State = np.concatenate((np.array([1,0]),State),axis=0)
+	
+	while diff > epsilon:
+		for i in range(0, nodes_num):
+			m = 2*i
+			i2 = i + nodes_num
+			type = system.nodes[i].type 
+			if 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:
+				z[m] = (np.real(system.nodes[i].power)*V[i].real + np.imag(system.nodes[i].power)*V[i].imag)/(np.abs(V[i])**2)
+				z[m+1] = (np.real(system.nodes[i].power)*V[i].imag - np.imag(system.nodes[i].power)*V[i].real)/(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:
+				z[m] = (np.real(system.nodes[i].power)*V[i].real + np.imag(system.nodes[i].power)*V[i].imag)/(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)
+		State = State + Delta_State
+		diff = np.amax(np.absolute(Delta_State))
+		
+		V = State[:nodes_num] + 1j * State[nodes_num:]
+				
+		num_iter = num_iter+1
+		
+	return V, num_iter
+	
+def calculateI(system, V):
+	"""
+	To calculate the branch currents
+	"""
+	Ymatrix, Adj = Ymatrix_calc(system)
+	I = np.zeros((len(system.branches)), dtype=np.complex)
+	for idx in range(len(system.branches)):
+		fr = system.branches[idx].start_node.index
+		to = system.branches[idx].end_node.index
+		I[idx] = - (V[fr] - V[to])*Ymatrix[fr][to]
+	
+	return I
+		
+def calculateInj(system, I):
+	"""
+	To calculate current injections at a node
+	"""
+	Iinj = np.zeros((len(system.nodes)), dtype=np.complex)
+	for k in range(0, (len(system.nodes))):
+		to=[]
+		fr=[]
+		for m in range(len(system.branches)):
+			if k==system.branches[m].start_node.index:
+				fr.append(m)
+			if k==system.branches[m].end_node.index:
+				to.append(m)
+		Iinj[k] = np.sum(I[to]) - np.sum(I[fr])
+	
+	return Iinj
+
+def calculateS1(system, V, I):
+	"""
+	To calculate powerflow on branches
+	"""
+	S1 = np.zeros((len(system.branches)), dtype=np.complex)
+	for i in range(0, len(system.branches)):
+		S1[i] = V[system.branches[i].start_node.index]*(np.conj(I[i]))
+	return S1
+	
+def calculateS2(system, V, I):
+	"""
+	To calculate powerflow on branches
+	"""
+	S2 = np.zeros((len(system.branches)), dtype=np.complex)
+	for i in range(0, len(system.branches)):
+		S2[i] = -V[system.branches[i].end_node.index]*(np.conj(I[i]))
+	return S2
+	
+def calculateSinj(V, Iinj):
+	"""
+	To calculate power injection at a node
+	"""
+	Sinj = np.multiply(V, np.conj(Iinj))
+	return Sinj
\ No newline at end of file

examples/CIGREMV_tests/comparison_dpsim_state_estimator.py

@@ -5,14 +5,15 @@ import copy
 import logging
 import matplotlib.pyplot as plt
 
-sys.path.append("../../cimpy")
+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
+sys.path.append("../../../dataprocessing")
+from villas.dataprocessing.readtools import *
 import nv_state_estimator_cim
 
 logging.basicConfig(filename='CIGRE.log', level=logging.INFO)
@@ -33,12 +34,13 @@ system = network.System()
 system.load_cim_data(res)
 Ymatr, Adj = network.Ymatrix_calc(system)
 
-loadflow_results = read_data(loadflow_results_file)
+#read Input-Ergebnisdatei
+ts_dpsim = read_timeseries_csv(loadflow_results_file)
 
 #order readed data according to system.nodes
-Vtrue=np.zeros(len(loadflow_results), dtype=np.complex_)
+Vtrue=np.zeros(len(ts_dpsim), dtype=np.complex_)
 for elem in range(len(system.nodes)): 
-	Vtrue[elem] = loadflow_results[system.nodes[elem].uuid]
+	Vtrue[elem] = ts_dpsim[system.nodes[elem].uuid].values[0]
 
 """ Write here the indexes of the nodes/branches where the measurements are"""
 V_idx = np.array([])

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