added new file nv_powerflow_cim.py and some examples

acs/state_estimation/network.py

 import numpy as np
+from enum import Enum
+
+class BusType(Enum):
+	SLACK = 1
+	slack = 1
+	PV = 2
+	pv = 2
+	PQ = 3
+	pq = 3
+
 
 class Node():
-	def __init__(self, name, uuid, v_mag, v_phase, p, q, index):
+	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
 		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.type = BusType[bus_type]
 				
 class Branch():
 	def __init__(self, r, x, start_node, end_node):
@@ -14,6 +25,8 @@ class Branch():
 		self.x = x
 		self.start_node = start_node
 		self.end_node = end_node
+		self.z = self.r + 1j*self.x
+		self.y = 1/self.z
 
 class System():	
 	def __init__(self):
@@ -31,7 +44,7 @@ class System():
 				self.P.append(value.pInjection)
 				self.Q.append(value.qInjection)
 				index=len(self.P)-1
-				self.node.append(Node(value.name, value.mRID, value.pInjection, value.qInjection, index))
+				self.nodes.append(Node(name=value.name, uuid=value.mRID, p=value.pInjection, q=value.qInjection, index=index))
 			elif value.__class__.__name__=="ACLineSegment":
 				length=value.length
 				if length==0.0:
@@ -40,25 +53,25 @@ class System():
 				self.bX.append(value.x*length)
 				startNode=res[value.startNodeID]
 				endNode=res[value.endNodeID]
-				self.branch.append(Branch(value.r, value.x, startNode, endNode))	
+				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)
 				startNode=res[value.startNodeID]
 				endNode=res[value.endNodeID]
-				self.branch.append(Branch(value.primaryConnection.r, value.primaryConnection.x, startNode, endNode))
+				self.branches.append(Branch(value.primaryConnection.r, value.primaryConnection.x, startNode, endNode))
 			else:
 				continue
 
 
-def load_python_data(nodes, branches):
+def load_python_data(nodes, branches, type):
 	system = System()
 	
 	for node_idx in range(0, nodes.num):
-		if node_idx == 0:			
-			system.nodes.append(Node('', '', nodes.P2[0], nodes.Q2[0], 0, 0, node_idx))
+		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:
-			system.nodes.append(Node('', '', 0, 0, nodes.P2[node_idx], nodes.Q2[node_idx], node_idx))
+			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]))
 	
 	for branch_idx in range(0, branches.num):
 		system.branches.append(Branch(branches.R[branch_idx], branches.X[branch_idx], 

acs/state_estimation/nv_powerflow_cim.py

+import numpy as np
+import math
+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
+		
+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
+		
+	Irx = np.zeros((branches_num),dtype=np.complex)
+	for idx in range(branches_num):
+		fr = system.branches[idx].start_node.index
+		to = system.branches[idx].end_node.index
+		Irx[idx] = - (V[fr] - V[to])*Ymatrix[fr][to]
+	Ir = np.real(Irx)
+	Ix = np.imag(Irx)
+	
+	I = Ir + 1j*Ix
+	Iinj_r = np.zeros(nodes_num)
+	Iinj_x = np.zeros(nodes_num)
+	
+	for k in range(0, nodes_num):
+		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_r[k] = np.sum(I[to].real) - np.sum(I[fr].real)
+		Iinj_x[k] = np.sum(I[to].imag) - np.sum(I[fr].imag)
+		
+	Iinj = Iinj_r + 1j * Iinj_x
+	Sinj_rx = np.multiply(V, np.conj(Iinj))
+	Sinj = np.real(Sinj_rx) + 1j * np.imag(Sinj_rx)
+	
+	S1=np.array([])
+	S2=np.array([])
+	I_conj=np.conj(I)
+	for i in range(0, branches_num):
+		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

examples/test_nv_powerflow.py

+import math
+import sys
+
+sys.path.append(r"C:\Users\Martin\Desktop\git\state-estimation\acs\state_estimation")
+sys.path.append(r"C:\Users\Martin\Desktop\git\state-estimation\examples")
+
+
+from network import *
+import py_95bus_network_data
+from nv_powerflow import *
+
+class PerUnit:
+    def __init__(self, S, V):
+        self.S = S
+        self.V = V
+        self.I = S/V
+        self.Z = S/(V**2)
+		
+		
+""" 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)
+V, I, Iinj, S1, S2, Sinj, num_iter = solve(branch, node)
+
+print("S1:")
+print(S1)
+print("\n\n")
+
+print("Sinj:")
+print(Sinj)
+print("\n\n")
+
+print("num_iter:" + str(num_iter))
+print("\n\n")
+
+

examples/test_nv_powerflow_cim.py

+import math
+import sys
+
+sys.path.append(r"C:\Users\Martin\Desktop\git\state-estimation\acs\state_estimation")
+sys.path.append(r"C:\Users\Martin\Desktop\git\state-estimation\examples")
+
+
+from network import *
+import py_95bus_network_data
+from nv_powerflow_cim import Ymatrix_calc as Ymatrix_calc_cim
+from nv_powerflow_cim import solve as solve_cim
+from nv_powerflow import *
+
+class PerUnit:
+    def __init__(self, S, V):
+        self.S = S
+        self.V = V
+        self.I = S/V
+        self.Z = S/(V**2)
+		
+		
+""" 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)
+system=load_python_data(node, branch, node.type)
+
+V, I, Iinj, S1, S2, Sinj, num_iter = solve_cim(system)
+
+print("S1_cim:")
+print(S1)
+print("\n\n")
+
+print("Sinj_cim:")
+print(Sinj)
+print("\n\n")
+
+print("num_iter:" + str(num_iter))
+print("\n\n")
+