Commit 707f91e4 authored by Jan Dinkelbach's avatar Jan Dinkelbach
Browse files

reformat code

parent 9b826343
This diff is collapsed.
import numpy as np
from enum import Enum
class BusType(Enum):
SLACK = 1
slack = 1
PV = 2
pv = 2
PQ = 3
pq = 3
SLACK = 1
slack = 1
PV = 2
pv = 2
PQ = 3
pq = 3
class Node():
def __init__(self, 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', pu=False):
self.uuid = uuid
self.index = index
self.baseVoltage = base_voltage
self.base_apparent_power = base_apparent_power
self.base_current = self.base_apparent_power/self.baseVoltage/np.sqrt(3)
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
def __init__(self, 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', pu=False):
self.uuid = uuid
self.index = index
self.baseVoltage = base_voltage
self.base_apparent_power = base_apparent_power
self.base_current = self.base_apparent_power / self.baseVoltage / np.sqrt(3)
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, uuid='', r=0.0, x=0.0, start_node=None, end_node=None,
base_voltage=1.0, base_apparent_power=1.0):
self.uuid = uuid
self.baseVoltage = base_voltage
self.base_apparent_power = base_apparent_power
self.base_current = self.base_apparent_power/self.baseVoltage/np.sqrt(3)
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")
def __init__(self, uuid='', r=0.0, x=0.0, start_node=None, end_node=None,
base_voltage=1.0, base_apparent_power=1.0):
self.uuid = uuid
self.baseVoltage = base_voltage
self.base_apparent_power = base_apparent_power
self.base_current = self.base_apparent_power / self.baseVoltage / np.sqrt(3)
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):
self.nodes = []
self.branches = []
self.Ymatrix = np.zeros([], dtype=np.complex)
self.Adjacencies = np.array([])
class System():
def __init__(self):
self.nodes=[]
self.branches=[]
self.Ymatrix = np.zeros([], dtype=np.complex)
self.Adjacencies = np.array([])
def load_cim_data(self, res, base_apparent_power):
"""
To fill the vectors node, branch, bR, bX, P and Q
"""
index = 0
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
node_type = self._getNodeType(element)
base_voltage = element.BaseVoltage.nominalVoltage
self.nodes.append(Node(uuid=uuid, base_voltage=base_voltage, v_mag=vmag,
base_apparent_power=base_apparent_power, v_phase=vphase,
p=pInj, q=qinj, index=index, bus_type=node_type))
index = index + 1
def load_cim_data(self, res, base_apparent_power):
"""
To fill the vectors node, branch, bR, bX, P and Q
"""
index = 0
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
node_type = self._getNodeType(element)
base_voltage = element.BaseVoltage.nominalVoltage
self.nodes.append(Node(uuid=uuid, base_voltage=base_voltage, v_mag=vmag,
base_apparent_power=base_apparent_power, v_phase=vphase,
p=pInj, q=qinj, index=index, bus_type=node_type))
index = index+1
for uuid, element in res.items():
if element.__class__.__name__=="ACLineSegment":
for node in self.nodes:
if element.startNodeID==node.uuid:
startNode = node
break
for node in self.nodes:
if element.endNodeID==node.uuid:
endNode=node
break
base_voltage = element.BaseVoltage.nominalVoltage
self.branches.append(Branch(uuid=uuid, r=element.r, x=element.x, start_node=startNode,
end_node=endNode, base_voltage=base_voltage,
base_apparent_power=base_apparent_power))
continue
elif element.__class__.__name__=="PowerTransformer":
for i in range(len(self.nodes)):
if element.startNodeID==self.nodes[i].uuid:
startNode=self.nodes[i]
break
for i in range(len(self.nodes)):
if element.endNodeID==self.nodes[i].uuid:
endNode=self.nodes[i]
break
#base voltage = high voltage side (=primaryConnection)
base_voltage = element.primaryConnection.BaseVoltage.nominalVoltage
self.branches.append(Branch(uuid=uuid, r=element.primaryConnection.r, x=element.primaryConnection.x,
start_node=startNode, end_node=endNode, base_voltage=base_voltage,
base_apparent_power=base_apparent_power))
continue
else:
continue
for uuid, element in res.items():
if element.__class__.__name__ == "ACLineSegment":
for node in self.nodes:
if element.startNodeID == node.uuid:
startNode = node
break
for node in self.nodes:
if element.endNodeID == node.uuid:
endNode = node
break
base_voltage = element.BaseVoltage.nominalVoltage
self.branches.append(Branch(uuid=uuid, r=element.r, x=element.x, start_node=startNode,
end_node=endNode, base_voltage=base_voltage,
base_apparent_power=base_apparent_power))
continue
elif element.__class__.__name__ == "PowerTransformer":
for i in range(len(self.nodes)):
if element.startNodeID == self.nodes[i].uuid:
startNode = self.nodes[i]
break
for i in range(len(self.nodes)):
if element.endNodeID == self.nodes[i].uuid:
endNode = self.nodes[i]
break
# base voltage = high voltage side (=primaryConnection)
base_voltage = element.primaryConnection.BaseVoltage.nominalVoltage
self.branches.append(Branch(uuid=uuid, r=element.primaryConnection.r, x=element.primaryConnection.x,
start_node=startNode, end_node=endNode, base_voltage=base_voltage,
base_apparent_power=base_apparent_power))
continue
else:
continue
#calculate admitance matrix
self.Ymatrix_calc()
# calculate admitance matrix
self.Ymatrix_calc()
def _getNodeType(self, node):
"""
return the type of a node: PQ, PV or SLACK
def _getNodeType(self, node):
"""
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"
@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.Adjacencies = [[] for _ in range(nodes_num)]
for branch in self.branches:
fr = branch.start_node.index
to = branch.end_node.index
self.Ymatrix[fr][to] -= branch.y_pu
self.Ymatrix[to][fr] -= branch.y_pu
self.Ymatrix[fr][fr] += branch.y_pu
self.Ymatrix[to][to] += branch.y_pu
self.Adjacencies[fr].append(to+1) #to + 1???
self.Adjacencies[to].append(fr+1) #fr + 1???
\ No newline at end of file
def Ymatrix_calc(self):
nodes_num = len(self.nodes)
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
to = branch.end_node.index
self.Ymatrix[fr][to] -= branch.y_pu
self.Ymatrix[to][fr] -= branch.y_pu
self.Ymatrix[fr][fr] += branch.y_pu
self.Ymatrix[to][to] += branch.y_pu
self.Adjacencies[fr].append(to + 1) # to + 1???
self.Adjacencies[to].append(fr + 1) # fr + 1???
......@@ -2,8 +2,9 @@ import numpy as np
from .network import BusType
from .results import Results
def solve(system):
"""It performs powerflow by using rectangular node voltage state variables and considering the current mismatch function.
"""It performs powerflow by using rectangular node voltage state variables and considering the current mismatch function.
Solve the non-linear powerflow problem stated by
......@@ -21,82 +22,85 @@ def solve(system):
H: Jacobian matrix
V: same as state but with complex numbers (i.e. [V0_re+j*V0_im, V1_re+j*V1_im, ...])
"""
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
node_type = system.nodes[i].type
if node_type == BusType.SLACK:
z[m] = np.real(system.nodes[i].voltage_pu)
z[m+1] = np.imag(system.nodes[i].voltage_pu)
H[m][i] = 1
H[m+1][i2] = 1
elif node_type is BusType.PQ:
H[m][i] = - np.real(system.Ymatrix[i][i])
H[m][i2] = np.imag(system.Ymatrix[i][i])
H[m+1][i] = - np.imag(system.Ymatrix[i][i])
H[m+1][i2] = - np.real(system.Ymatrix[i][i])
idx1 = np.subtract(system.Adjacencies[i],1)
idx2 = idx1 + nodes_num
H[m][idx1] = - np.real(system.Ymatrix[i][idx1])
H[m][idx2] = np.imag(system.Ymatrix[i][idx1])
H[m+1][idx1] = - np.imag(system.Ymatrix[i][idx1])
H[m+1][idx2] = - np.real(system.Ymatrix[i][idx1])
elif node_type is BusType.PV:
z[m+1] = np.real(system.nodes[i].power)
H[m][i] = - np.real(system.Ymatrix[i][i])
H[m][i2] = np.imag(system.Ymatrix[i][i])
idx1 = np.subtract(system.Adjacencies[i],1)
idx2 = idx1 + nodes_num
H[m][idx1] = - np.real(system.Ymatrix[i][idx1])
H[m][idx2] = np.imag(system.Ymatrix[i][idx1])
epsilon = 10**(-10)
diff = 5
V = np.ones(nodes_num) + 1j* np.zeros(nodes_num)
num_iter = 0
state = np.concatenate((np.ones(nodes_num),np.zeros(nodes_num)),axis=0)
while diff > epsilon:
for i in range(0, nodes_num):
m = 2*i
i2 = i + nodes_num
node_type = system.nodes[i].type
if node_type is BusType.SLACK:
h[m] = np.inner(H[m],state)
h[m+1] = np.inner(H[m+1],state)
elif node_type is BusType.PQ:
z[m] = (np.real(system.nodes[i].power_pu)*np.real(V[i]) + np.imag(system.nodes[i].power_pu)*np.imag(V[i]))/(np.abs(V[i])**2)
z[m+1] = (np.real(system.nodes[i].power_pu)*np.imag(V[i]) - np.imag(system.nodes[i].power_pu)*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)
elif node_type is BusType.PV:
z[m] = (np.real(system.nodes[i].power_pu)*np.real(V[i])+ np.imag(system.nodes[i].power_pu)*np.imag(V[i]))(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
# calculate all the other quantities of the grid
powerflow_results = Results(system)
powerflow_results.load_voltages(V)
powerflow_results.calculate_all()
return powerflow_results, num_iter
\ No newline at end of file
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
node_type = system.nodes[i].type
if node_type == BusType.SLACK:
z[m] = np.real(system.nodes[i].voltage_pu)
z[m + 1] = np.imag(system.nodes[i].voltage_pu)
H[m][i] = 1
H[m + 1][i2] = 1
elif node_type is BusType.PQ:
H[m][i] = - np.real(system.Ymatrix[i][i])
H[m][i2] = np.imag(system.Ymatrix[i][i])
H[m + 1][i] = - np.imag(system.Ymatrix[i][i])
H[m + 1][i2] = - np.real(system.Ymatrix[i][i])
idx1 = np.subtract(system.Adjacencies[i], 1)
idx2 = idx1 + nodes_num
H[m][idx1] = - np.real(system.Ymatrix[i][idx1])
H[m][idx2] = np.imag(system.Ymatrix[i][idx1])
H[m + 1][idx1] = - np.imag(system.Ymatrix[i][idx1])
H[m + 1][idx2] = - np.real(system.Ymatrix[i][idx1])
elif node_type is BusType.PV:
z[m + 1] = np.real(system.nodes[i].power)
H[m][i] = - np.real(system.Ymatrix[i][i])
H[m][i2] = np.imag(system.Ymatrix[i][i])
idx1 = np.subtract(system.Adjacencies[i], 1)
idx2 = idx1 + nodes_num
H[m][idx1] = - np.real(system.Ymatrix[i][idx1])
H[m][idx2] = np.imag(system.Ymatrix[i][idx1])
epsilon = 10 ** (-10)
diff = 5
V = np.ones(nodes_num) + 1j * np.zeros(nodes_num)
num_iter = 0
state = np.concatenate((np.ones(nodes_num), np.zeros(nodes_num)), axis=0)
while diff > epsilon:
for i in range(0, nodes_num):
m = 2 * i
i2 = i + nodes_num
node_type = system.nodes[i].type
if node_type is BusType.SLACK:
h[m] = np.inner(H[m], state)
h[m + 1] = np.inner(H[m + 1], state)
elif node_type is BusType.PQ:
z[m] = (np.real(system.nodes[i].power_pu) * np.real(V[i]) + np.imag(system.nodes[i].power_pu) * np.imag(
V[i])) / (np.abs(V[i]) ** 2)
z[m + 1] = (np.real(system.nodes[i].power_pu) * np.imag(V[i]) - np.imag(
system.nodes[i].power_pu) * 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)
elif node_type is BusType.PV:
z[m] = (np.real(system.nodes[i].power_pu) * np.real(V[i]) + np.imag(system.nodes[i].power_pu) * np.imag(
V[i]))(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
# calculate all the other quantities of the grid
powerflow_results = Results(system)
powerflow_results.load_voltages(V)
powerflow_results.calculate_all()
return powerflow_results, num_iter
This diff is collapsed.
......@@ -2,248 +2,250 @@ import cmath
import numpy as np
from villas.dataprocessing.readtools import read_timeseries_dpsim
class ResultsNode():
def __init__(self, topo_node):
self.topology_node = topo_node
self.voltage = complex(0, 0)
self.current = complex(0, 0)
self.power = complex(0, 0)
self.voltage_pu = complex(0, 0)
self.current_pu = complex(0, 0)
self.power_pu = complex(0, 0)
def __init__(self, topo_node):
self.topology_node = topo_node
self.voltage = complex(0, 0)
self.current = complex(0, 0)
self.power = complex(0, 0)
self.voltage_pu = complex(0, 0)
self.current_pu = complex(0, 0)
self.power_pu = complex(0, 0)
class ResultsBranch():
def __init__(self, topo_branch):
self.topology_branch = topo_branch
self.current = complex(0, 0)
self.power = complex(0, 0) #complex power flow at branch, measured at initial node
self.power2 = complex(0, 0) #complex power flow at branch, measured at final node
self.current_pu = complex(0, 0)
self.power_pu = complex(0, 0) #complex power flow at branch, measured at initial node
self.power2_pu = complex(0, 0) #complex power flow at branch, measured at final node
class Results():
def __init__(self, system):
self.nodes=[]
self.branches=[]
self.Ymatrix=system.Ymatrix
self.Adjacencies=system.Adjacencies
for node in system.nodes:
self.nodes.append(ResultsNode(topo_node=node))
for branch in system.branches:
self.branches.append(ResultsBranch(topo_branch=branch))
def read_data_dpsim(self, file_name, pu=False):
"""
read the voltages from a dpsim input file
@param file_name
@param pu: - True if voltages are expressed in per unit system
"""
loadflow_results = read_timeseries_dpsim(file_name, print_status=False)
if pu==True:
for node in self.nodes:
node.voltage_pu = loadflow_results[node.topology_node.uuid].values[0]
node.voltage = node.voltage_pu*node.topology_node.baseVoltage
elif pu==False:
for node in self.nodes:
node.voltage = loadflow_results[node.topology_node.uuid].values[0]
node.voltage_pu = node.voltage/node.topology_node.baseVoltage
def load_voltages(self, V):
"""
load the voltages of V-array (result of powerflow_cim.solve)
"""
for index in range(len(V)):
for node in self.nodes:
if node.topology_node.index == index:
node.voltage_pu = V[index]
node.voltage = node.voltage_pu*node.topology_node.baseVoltage
def calculate_all(self):
"""
calculate all quantities of the grid
"""
self.calculateI()
self.calculateIinj()
self.calculateSinj()
self.calculateIinj()
self.calculateS1()
self.calculateS2()
def calculateI(self):
"""
To calculate the branch currents
"""
for branch in self.branches:
fr = branch.topology_branch.start_node.index
to = branch.topology_branch.end_node.index
branch.current_pu = - (self.nodes[fr].voltage_pu - self.nodes[to].voltage_pu)*self.Ymatrix[fr][to]
branch.current = branch.current_pu*branch.topology_branch.base_current
def calculateIinj(self):
"""
calculate current injections at a node
"""
for node in self.nodes:
to=complex(0, 0) #sum of the currents flowing to the node
fr=complex(0, 0) #sum of the currents flowing from the node
for branch in self.branches:
if node.topology_node.index==branch.topology_branch.start_node.index:
fr = fr + branch.current_pu
if node.topology_node.index==branch.topology_branch.end_node.index:
to = to + branch.current_pu
node.current_pu = to - fr
node.current = node.current_pu*node.topology_node.base_current
def calculateSinj(self):
"""
calculate power injection at a node