Aufgrund einer Störung des s3 Storage, könnten in nächster Zeit folgende GitLab Funktionen nicht zur Verfügung stehen: LFS, Container Registry, Job Artifacs, Uploads (Wiki, Bilder, Projekt-Exporte). Wir bitten um Verständnis. Es wird mit Hochdruck an der Behebung des Problems gearbeitet. Weitere Informationen zur Störung des Object Storage finden Sie hier: https://maintenance.itc.rwth-aachen.de/ticket/status/messages/59-object-storage-pilot

Aufgrund einer Wartung wird GitLab am 03.08. zwischen 8:00 und 9:00 Uhr kurzzeitig nicht zur Verfügung stehen. / Due to maintenance, GitLab will be temporarily unavailable on 03.08. between 8:00 and 9:00 am.

Commit b6070eb6 authored by Jan Dinkelbach's avatar Jan Dinkelbach
Browse files

Merge branch 'new-cimpy'

parents 4939ef5b 9d909403
import logging
import numpy as np
from enum import Enum
class BusType(Enum):
SLACK = 1
slack = 1
......@@ -12,18 +12,20 @@ class BusType(Enum):
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):
def __init__(self, uuid='', name='', 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, ideal_connected_with=''):
self.uuid = uuid
self.name = name
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["PQ"]
self.voltage = v_mag * np.cos(np.radians(v_phase)) + 1j * v_mag * np.sin(np.radians(v_phase))
self.voltage = complex(v_mag * np.cos(np.radians(v_phase)), v_mag * np.sin(np.radians(v_phase)))
self.power = complex(p, q)
self.power_pu = complex(p, q) / self.base_apparent_power
self.voltage_pu = self.voltage / self.baseVoltage
self.ideal_connected_with = ideal_connected_with
def __str__(self):
str = 'class=Node\n'
......@@ -31,6 +33,7 @@ class Node():
for key in attributes.keys():
str = str + key + '={}\n'.format(attributes[key])
return str
class Branch():
def __init__(self, uuid='', r=0.0, x=0.0, start_node=None, end_node=None,
......@@ -58,121 +61,221 @@ class Branch():
str = str + key + '={}\n'.format(attributes[key])
return str
class Breaker():
def __init__(self, from_node, to_node, is_open=True):
"""
:param from_node:
:param to_node:
:param is_open: True if the breaker is considered open and False if the broker is closed
"""
self.from_node = from_node
self.to_node = to_node
self.is_open = is_open
def __str__(self):
str = 'class=Breaker\n'
attributes = self.__dict__
for key in attributes.keys():
str = str + key + '={}\n'.format(attributes[key])
return str
def open_breaker(self):
self.is_open == True
self.to_node.ideal_connected_with = ''
def close_breaker(self):
self.is_open == False
self.to_node.ideal_connected_with = self.from_node.uuid
class System():
def __init__(self):
self.nodes = []
self.branches = []
self.breakers = []
self.Ymatrix = np.zeros([], dtype=np.complex)
self.Adjacencies = np.array([])
def get_node_by_uuid(self, node_uuid):
for node in self.nodes:
if node.uuid == node_uuid:
return node
return False
def get_node_by_index(self, index):
"""
return the node with node.index==index
"""
for node in self.nodes:
if (node.index==index) and (node.ideal_connected_with=='') :
return node
return None
def get_nodes_num(self):
"""
return the number of nodes in the list system.nodes
Warning: if any node is ideally connected to another node,
the counter is increased only one time
"""
nodes_num=0
for node in self.nodes:
if node.ideal_connected_with=='':
nodes_num+=1
return nodes_num
def reindex_nodes_list(self):
"""
Reenumerate the nodes in system.nodes
If any node is ideally connected to another node,
both receive the same index
"""
index = 0
remaining_nodes_list = []
for node in self.nodes:
if node.ideal_connected_with=='':
node.index=index
index+=1
else:
remaining_nodes_list.append(node)
for node in remaining_nodes_list:
node.index = self.get_node_by_uuid(node.ideal_connected_with).index
def load_cim_data(self, res, base_apparent_power):
"""
To fill the vectors node, branch, bR, bX, P and Q
fill the vectors node, branch and breakers
"""
index = 0
list_TPNode = [elem for elem in res.values() if elem.__class__.__name__ == "TopologicalNode"]
list_SvVoltage = [elem for elem in res.values() if elem.__class__.__name__ == "SvVoltage"]
list_SvPowerFlow = [elem for elem in res.values() if elem.__class__.__name__ == "SvPowerFlow"]
list_EnergySources = [elem for elem in res.values() if elem.__class__.__name__ == "EnergySource"]
list_EnergyConsumer = [elem for elem in res.values() if elem.__class__.__name__ == "EnergyConsumer"]
list_ACLineSegment = [elem for elem in res.values() if elem.__class__.__name__ == "ACLineSegment"]
list_PowerTransformer = [elem for elem in res.values() if elem.__class__.__name__ == "PowerTransformer"]
list_Terminals = [elem for elem in res.values() if elem.__class__.__name__ == "Terminal"]
list_Terminals_ES = [elem for elem in list_Terminals if elem.ConductingEquipment.__class__.__name__ == "EnergySource"]
list_Terminals_EC = [elem for elem in list_Terminals if elem.ConductingEquipment.__class__.__name__ == "EnergyConsumer"]
list_PowerTransformerEnds = [elem for elem in res.values() if elem.__class__.__name__ == "PowerTransformerEnd"]
list_Breakers = [elem for elem in res.values() if elem.__class__.__name__ == "Breaker"]
#create nodes
for TPNode in list_TPNode:
uuid_TPNode = TPNode.mRID
name = TPNode.name
vmag = 0.0
vphase = 0.0
pInj = 0.0
qinj = 0.0
qInj = 0.0
for obj_SvVoltage in list_SvVoltage:
if obj_SvVoltage.TopologicalNode[0].mRID == uuid_TPNode:
if obj_SvVoltage.TopologicalNode.mRID == uuid_TPNode:
vmag = obj_SvVoltage.v
vphase = obj_SvVoltage.angle
break
for obj_SvPowerFlow in list_SvPowerFlow:
if obj_SvPowerFlow.Terminal[0].TopologicalNode[0].mRID == uuid_TPNode:
if obj_SvPowerFlow.Terminal.TopologicalNode.mRID == uuid_TPNode:
pInj += obj_SvPowerFlow.p
qinj += obj_SvPowerFlow.q
base_voltage = TPNode.BaseVoltage[0].nominalVoltage
self.nodes.append(Node(uuid=uuid_TPNode, base_voltage=base_voltage, v_mag=vmag,
qInj += obj_SvPowerFlow.q
for obj_Terminal in list_Terminals_ES:
if obj_Terminal.TopologicalNode.mRID == uuid_TPNode:
for obj_EnergySource in list_EnergySources:
if obj_EnergySource.mRID == obj_Terminal.ConductingEquipment.mRID:
pInj += obj_EnergySource.activePower
qInj += obj_EnergySource.reactivePower
for obj_Terminal in list_Terminals_EC:
if obj_Terminal.TopologicalNode.mRID == uuid_TPNode:
for obj_EnergyConsumer in list_EnergyConsumer:
if obj_EnergyConsumer.mRID == obj_Terminal.ConductingEquipment.mRID:
pInj += obj_EnergyConsumer.p
qInj += obj_EnergyConsumer.q
base_voltage = TPNode.BaseVoltage.nominalVoltage
self.nodes.append(Node(name=name, uuid=uuid_TPNode, base_voltage=base_voltage, v_mag=vmag,
base_apparent_power=base_apparent_power, v_phase=vphase,
p=pInj, q=qinj, index=index))
p=pInj, q=qInj, index=index))
index = index + 1
#create branches ACLineSegment
self._setNodeType(list_Terminals)
#create branches type ACLineSegment
for ACLineSegment in list_ACLineSegment:
uuid_ACLineSegment = ACLineSegment.mRID
start_node_id, end_node_id = self._get_node_uuids(list_Terminals, uuid_ACLineSegment)
start_node = None
end_node = None
nodes = self._get_nodes(list_Terminals, uuid_ACLineSegment)
start_node = nodes[0]
end_node = nodes[1]
for node in self.nodes:
if start_node_id == node.uuid:
start_node = node
elif end_node_id == node.uuid:
end_node = node
base_voltage = ACLineSegment.BaseVoltage[0].nominalVoltage
base_voltage = ACLineSegment.BaseVoltage.nominalVoltage
self.branches.append(Branch(uuid=uuid_ACLineSegment, r=ACLineSegment.r, x=ACLineSegment.x,
start_node=start_node, end_node=end_node,
base_voltage=base_voltage, base_apparent_power=base_apparent_power))
#create branches power transformer
#create branches type powerTransformer
for power_transformer in list_PowerTransformer:
uuid_power_transformer = power_transformer.mRID
start_node_id, end_node_id = self._get_node_uuids(list_Terminals, uuid_power_transformer)
start_node = None
end_node = None
for node in self.nodes:
if start_node_id == node.uuid:
start_node = node
elif end_node_id == node.uuid:
end_node = node
nodes = self._get_nodes(list_Terminals, uuid_power_transformer)
start_node = nodes[0]
end_node = nodes[1]
# base voltage = high voltage side (=primaryConnection)
primary_connection = self._get_primary_connection(list_PowerTransformerEnds, uuid_power_transformer)
base_voltage = primary_connection.BaseVoltage[0].nominalVoltage
base_voltage = primary_connection.BaseVoltage.nominalVoltage
self.branches.append(Branch(uuid=uuid_power_transformer, r=primary_connection.r, x=primary_connection.x,
start_node=start_node, end_node=end_node, base_voltage=base_voltage,
base_apparent_power=base_apparent_power))
self._setNodeType(list_Terminals)
# calculate admitance matrix
#create breakers
for obj_Breaker in list_Breakers:
is_open = obj_Breaker.normalOpen
nodes = self._get_nodes(list_Terminals, obj_Breaker.mRID)
self.breakers.append(Breaker(from_node=nodes[0], to_node=nodes[1], is_open=is_open))
#if the breaker is open == closed --> close broker
if is_open == False:
self.breakers[-1].close_breaker(self)
else:
self.breakers[-1].ideal_connected_with = ''
#calculate admitance matrix
self.Ymatrix_calc()
def _get_node_uuids(self, list_Terminals, elem_uuid):
def _get_nodes(self, list_Terminals, elem_uuid):
"""
get the uuids of the start and end nodes of the element with uuid=elem_uuid
This function can only used with the elements ACLineSegment, PowerTransformer, Switch and EnergyConsumer
get the the start and end node of the element with uuid=elem_uuid
This function can used only with element which are connected
to 2 topologicalNodes, for example: ACLineSegment, PowerTransformer and Breaker
:param list_Terminals: list of all elements of type Terminal
:param elem_uuid: uuid of the element for which the start and end node ID are searched
:return tuple: (startNodeID, endNodeID)
:return list: [startNodeID, endNodeID]
"""
start_node_uuid = ""
end_node_uuid = ""
list_elements = ["ACLineSegment", "PowerTransformer", "Switch", "EnergyConsumer"]
start_node_uuid = None
end_node_uuid = None
for terminal in list_Terminals:
#TODO: has always the list ConductingEquipment only one element?????
conductingEquipment = terminal.ConductingEquipment[0]
if (conductingEquipment.mRID != elem_uuid):
if (terminal.ConductingEquipment.mRID != elem_uuid):
continue
elif conductingEquipment.__class__.__name__ in list_elements:
sequence_number = terminal.sequenceNumber
#TODO: has always the list Terminal.TopologicalNode only one element?????
topological_node = terminal.TopologicalNode[0]
if sequence_number == 1:
start_node_uuid = topological_node.mRID
elif sequence_number == 2:
end_node_uuid = topological_node.name
sequence_number = terminal.sequenceNumber
if sequence_number == 1:
start_node_uuid = terminal.TopologicalNode.mRID
elif sequence_number == 2:
end_node_uuid = terminal.TopologicalNode.mRID
return (start_node_uuid, end_node_uuid)
start_node = None
end_node = None
if (start_node_uuid==None):
print('WARNING: It could not find a start node for the element with uuid={}'.format(elem_uuid))
else:
start_node = self.get_node_by_uuid(start_node_uuid)
if (end_node_uuid==None):
print('WARNING: It could not find a end node for the element with uuid={}'.format(elem_uuid))
else:
end_node = self.get_node_by_uuid(end_node_uuid)
return [start_node, end_node]
def _get_primary_connection(self, list_PowerTransformerEnds, elem_uuid):
"""
......@@ -180,49 +283,57 @@ class System():
:param list_PowerTransformerEnd: list of all elements of type PowerTransformerEnd
:param elem_uuid: uuid of the power transformer for which the primary connection is searched
:return: primary_connection
"""
"""
primary_connection = None
power_transformer_ends = []
#search for two elements of class powertransformerend that point to the powertransformer with ID = elem_uuid
for power_transformer_end in list_PowerTransformerEnds:
if power_transformer_end.PowerTransformer[0].mRID == elem_uuid:
power_transformer = None
if isinstance(power_transformer_end.PowerTransformer, list):
if (len(power_transformer_end.PowerTransformer)!=1):
print('WARNING: len(power_transformer_end.PowerTransformer)!=1 for the element with uuid={}. \
The first element will be used'.format(power_transformer_end.mRID))
power_transformer = power_transformer_end.PowerTransformer[0]
else:
power_transformer = power_transformer_end.PowerTransformer
if power_transformer.mRID == elem_uuid:
power_transformer_ends.append(power_transformer_end)
if power_transformer_ends[0].BaseVoltage[0].nominalVoltage>=power_transformer_ends[1].BaseVoltage[0].nominalVoltage:
if power_transformer_ends[0].BaseVoltage.nominalVoltage>=power_transformer_ends[1].BaseVoltage.nominalVoltage:
primary_connection=power_transformer_ends[0]
elif power_transformer_ends[1].BaseVoltage[0].nominalVoltage>=power_transformer_ends[0].BaseVoltage[0].nominalVoltage:
elif power_transformer_ends[1].BaseVoltage.nominalVoltage>=power_transformer_ends[0].BaseVoltage.nominalVoltage:
primary_connection=power_transformer_ends[1]
return primary_connection
def _setNodeType(self, list_Terminals):
"""
set the parameter type of all elements of the list self.nodes
set the parameter "type" of all elements of the list self.nodes
:param list_PowerTransformerEnd: list of all elements of type Terminal
:return None
"""
#get a list of Terminals for which the ConductingEquipment is a element of class ExternalNetworkInjection
list_Terminals_ENI = [elem for elem in list_Terminals if elem.ConductingEquipment[0].__class__.__name__ == "ExternalNetworkInjection"]
list_Terminals_ENI = [elem for elem in list_Terminals if elem.ConductingEquipment.__class__.__name__ == "ExternalNetworkInjection"]
for terminal in list_Terminals_ENI:
#TODO: is it correct to use the element 0 for it? Is len(terminal.TopologicalNode) greater than one?
node_uuid = terminal.TopologicalNode[0].mRID
node_uuid = terminal.TopologicalNode.mRID
for node in self.nodes:
if node.uuid == node_uuid:
node.type = BusType["SLACK"]
#TODO the search for PV nodes has not been tested yet
#get a list of Terminals for which the ConductingEquipment is a element of class SynchronousMachine
list_Terminals_SM = [elem for elem in list_Terminals if elem.ConductingEquipment[0].__class__.__name__ == "SynchronousMachine"]
list_Terminals_SM = [elem for elem in list_Terminals if elem.ConductingEquipment.__class__.__name__ == "SynchronousMachine"]
for terminal in list_Terminals_SM:
#TODO: is it correct to use the element 0 for it? Is len(terminal.TopologicalNode) greater than one?
node_uuid = terminal.TopologicalNode[0].MRID
node_uuid = terminal.TopologicalNode.mRID
for node in self.nodes:
if node.uuid == node_uuid:
node.type = BusType["PV"]
def Ymatrix_calc(self):
nodes_num = len(self.nodes)
self.reindex_nodes_list()
nodes_num = self.get_nodes_num()
self.Ymatrix = np.zeros((nodes_num, nodes_num), dtype=np.complex)
self.Adjacencies = [[] for _ in range(nodes_num)]
for branch in self.branches:
......@@ -234,3 +345,17 @@ class System():
self.Ymatrix[to][to] += branch.y_pu
self.Adjacencies[fr].append(to + 1) # to + 1???
self.Adjacencies[to].append(fr + 1) # fr + 1???
#Testing functions
def print_nodes_names(self):
for node in self.nodes:
print('{} {}'.format(node.name, node.index))
def print_node_types(self):
for node in self.nodes:
print('{} {}'.format(node.name, node.type))
def print_power(self):
for node in self.nodes:
print('{} {}'.format(node.name, node.power))
......@@ -23,43 +23,44 @@ def solve(system):
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)
nodes_num = system.get_nodes_num()
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])
for node in system.nodes:
if node.ideal_connected_with =='':
i = node.index
m = 2 * i
i2 = i + nodes_num
node_type = node.type
if node_type == BusType.SLACK:
z[m] = np.real(node.voltage_pu)
z[m + 1] = np.imag(node.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(node.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)
#epsilon = 0.01
diff = 5
V = np.ones(nodes_num) + 1j * np.zeros(nodes_num)
num_iter = 0
......@@ -67,27 +68,29 @@ def solve(system):
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]))
for node in system.nodes:
if node.ideal_connected_with =='':
i = node.index
m = 2 * i
i2 = i + nodes_num
node_type = node.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(node.power_pu) * np.real(V[i]) +
np.imag(node.power_pu) * np.imag(V[i])) / (np.abs(V[i]) ** 2)
z[m + 1] = (np.real(node.power_pu) * np.imag(V[i]) -
np.imag(node.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(node.power_pu) * np.real(V[i]) +
np.imag(node.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)
......
......@@ -15,7 +15,8 @@ def DsseCall(system, measurements):
"""
# select type of Estimator.
# if at least a PMU is present we launch the combined estimator, otherwise a simple traditional etimator
# 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:
......@@ -29,7 +30,7 @@ def DsseCall(system, measurements):
est_code = trad_code + PMU_code
# number of nodes of the grid
nodes_num = len(system.nodes)
nodes_num = system.get_nodes_num()
Gmatrix = system.Ymatrix.real
Bmatrix = system.Ymatrix.imag
......@@ -39,7 +40,7 @@ def DsseCall(system, measurements):
# run Estimator.
if est_code == 1:
Vest = DsseTrad(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix)
Vest = DsseTrad(nodes_num, measurements, Gmatrix, Bmatrix, Yabs_matrix, Yphase_matrix, Adj)
elif est_code == 2:
Vest = DssePmu(nodes_num, measurements, Gmatrix, Bmatrix, Adj)
else:
......@@ -53,7 +54,7 @@ def DsseCall(system, measurements):
return results