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

Merge branch 'development'

parents f52e43df c6f7e7a2
......@@ -26,7 +26,7 @@ for file in xml_files:
xml_files_abs.append(os.path.abspath(file))
# read cim files and create new network.Systen object
res, _ = cimpy.cim_import(xml_files_abs, "cgmes_v2_4_15")
res, _, _ = cimpy.cim_import(xml_files_abs, "cgmes_v2_4_15")
system = network.System()
base_apparent_power = 25 # MW
system.load_cim_data(res, base_apparent_power)
......@@ -40,7 +40,7 @@ V_unc = 0
I_unc = 0
Sinj_unc = 0
S_unc = 0
Pmu_mag_unc = 1
Pmu_mag_unc = 0
Pmu_phase_unc = 0
# Create measurements data structures
......
......@@ -63,13 +63,23 @@ class MeasurementSet:
to update meas_value of a specific measurement object in the measurements array
"""
# only update measurements that are already included in the measurements set
for meas in self.measurements:
if meas.element.uuid == element_uuid and meas.meas_type == meas_type:
if not value_in_pu:
if meas.meas_type == MeasType.Vpmu_mag:
meas_value = meas_data / (meas.element.baseVoltage * 1000 / np.sqrt(
3)) # TODO - Fix phase-to-phase voltage problem
# special case for voltage magnitude: SOGNO interface only knows Vpmu_mag while measurement set distincts between Vpmu_mag and V_mag
if meas.element.uuid == element_uuid and meas_type == MeasType.Vpmu_mag and (meas.meas_type == MeasType.Vpmu_mag or meas.meas_type == MeasType.V_mag):
# volt pu conversion assuming that meas_data from device are in volts and single-phase value according to sogno interface
# while baseVoltage from CIM in [kV] and three-phase value
if not value_in_pu and (meas.meas_type == MeasType.Vpmu_mag or meas.meas_type == MeasType.V_mag):
meas_value = meas_data / (meas.element.baseVoltage / np.sqrt(3) * 1000)
meas.meas_value = meas_value
# case for other measurements
elif meas.element.uuid == element_uuid and meas.meas_type == meas_type:
# power pu conversion assuming that meas_data from device are in watts and single-phase value according to sogno interface
# while baseApparent power in [MW] and three-phase value
if not value_in_pu and (meas.meas_type == MeasType.S1_real or meas.meas_type == MeasType.S1_imag):
meas_value = meas_data / (meas.element.base_apparent_power / 3 * 1e6)
meas.meas_value = meas_value
def read_measurements_from_file(self, powerflow_results, file_name):
"""
......@@ -161,7 +171,7 @@ class MeasurementSet:
self.create_measurement(element, ElemType.Branch, MeasType.Ipmu_mag, meas_value_ideal_mag, unc_mag)
self.create_measurement(element, ElemType.Branch, MeasType.Ipmu_phase, meas_value_ideal_phase, unc_phase)
def meas_creation(self, dist="normal", seed=None):
def meas_creation(self, dist="normal", seed=None, type="simulation"):
"""
It calculates the measured values (affected by uncertainty) at the measurement points
which distribution should be used? if gaussian --> stddev must be divided by 3
......@@ -171,23 +181,27 @@ class MeasurementSet:
@param seed: - normal: use normal distribution (-->std_dev are divided by 3)
- uniform: use normal distribution
"""
np.random.seed(seed)
if dist == "normal":
err_pu = np.random.normal(0, 1, len(self.measurements))
for index, measurement in enumerate(self.measurements):
if measurement.meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.meas_value_ideal * measurement.std_dev
elif measurement.meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.std_dev
measurement.meas_value = measurement.meas_value_ideal + zdev * err_pu[index]
elif dist == "uniform":
err_pu = np.random.uniform(-1, 1, len(self.measurements))
for index, measurement in enumerate(self.measurements):
if measurement.meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = (measurement.meas_value_ideal * measurement.std_dev)
elif measurement.meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.std_dev
measurement.meas_value = measurement.meas_value_ideal + np.multiply(3 * zdev, err_pu[index])
if type == "simulation":
np.random.seed(seed)
if dist == "normal":
err_pu = np.random.normal(0, 0, len(self.measurements))
for index, measurement in enumerate(self.measurements):
if measurement.meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.meas_value_ideal * measurement.std_dev
elif measurement.meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.std_dev
measurement.meas_value = measurement.meas_value_ideal + zdev * err_pu[index]
elif dist == "uniform":
err_pu = np.random.uniform(-1, 1, len(self.measurements))
for index, measurement in enumerate(self.measurements):
if measurement.meas_type not in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = (measurement.meas_value_ideal * measurement.std_dev)
elif measurement.meas_type in [MeasType.Ipmu_phase, MeasType.Vpmu_phase]:
zdev = measurement.std_dev
measurement.meas_value = measurement.meas_value_ideal + np.multiply(3 * zdev, err_pu[index])
elif type == "field":
for measurement in self.measurements:
measurement.meas_value = measurement.meas_value_ideal
def meas_creation_test(self, err_pu):
"""
......@@ -204,7 +218,7 @@ class MeasurementSet:
# measurement.meas_value = measurement.meas_value_ideal + measurement.std_dev*err_pu[index]
def getMeasurements(self, type):
def getMeasurementsOfType(self, type):
"""
return an array with all measurements of type "type" in the array MeasurementSet.measurements.
"""
......@@ -283,6 +297,21 @@ class MeasurementSet:
meas_real[ipmu_phase_index] = iamp * np.sin(itheta)
return meas_real
def getSortedMeasurementSet(self):
"""
Sorts measurements in the order required by the SE algorithm
"""
sortedMeasurementSet = MeasurementSet()
# Sort measurements in the order required by the SE algorithm
# Required order: Vmag, Pinj, Qinj, P1, Q1, P2, Q2, Imag, Vpmu_mag, Vpmu_phase, Ipmu_mag, Ipmu_phase
for type_meas in [MeasType.V_mag, MeasType.Sinj_real, MeasType.Sinj_imag, MeasType.S1_real, MeasType.S1_imag, \
MeasType.S2_real, MeasType.S2_imag, MeasType.I_mag, MeasType.Vpmu_mag, MeasType.Vpmu_phase, MeasType.Ipmu_mag, MeasType.Ipmu_phase]:
sortedMeasurementSet.measurements += self.getMeasurementsOfType(type_meas)
return sortedMeasurementSet
def getStd_Dev(self):
"""
......
......@@ -95,7 +95,6 @@ class System():
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:
......@@ -179,8 +178,8 @@ class System():
break
for obj_SvPowerFlow in list_SvPowerFlow:
if obj_SvPowerFlow.Terminal.TopologicalNode.mRID == uuid_TPNode:
pInj += obj_SvPowerFlow.p
qInj += obj_SvPowerFlow.q
pInj -= obj_SvPowerFlow.p
qInj -= obj_SvPowerFlow.q
for obj_Terminal in list_Terminals_ES:
if obj_Terminal.TopologicalNode.mRID == uuid_TPNode:
for obj_EnergySource in list_EnergySources:
......@@ -191,8 +190,8 @@ class System():
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
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,
......@@ -335,16 +334,16 @@ class System():
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)]
self.Bmatrix = np.zeros((nodes_num, nodes_num), dtype=np.complex)
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???
self.Ymatrix[fr][fr] += branch.y_pu # + branch.b_pu
self.Ymatrix[to][to] += branch.y_pu # + branch.b_pu
# self.Bmatrix[fr][to] += branch.b_pu
# self.Bmatrix[to][fr] += branch.b_pu
#Testing functions
def print_nodes_names(self):
......
......@@ -40,24 +40,14 @@ def solve(system):
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])
H[m][:nodes_num] = np.real(system.Ymatrix[i])
H[m][nodes_num:] = - np.imag(system.Ymatrix[i])
H[m+1][:nodes_num] = np.imag(system.Ymatrix[i])
H[m+1][nodes_num:] = np.real(system.Ymatrix[i])
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])
H[m][:nodes_num] = np.real(system.Ymatrix[i])
H[m][nodes_num:] = - np.imag(system.Ymatrix[i])
epsilon = 10 ** (-10)
#epsilon = 0.01
......@@ -86,7 +76,7 @@ def solve(system):
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)
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]))
......
This diff is collapsed.
......@@ -44,7 +44,7 @@ class Results():
self.nodes = []
self.branches = []
self.Ymatrix = system.Ymatrix
self.Adjacencies = system.Adjacencies
self.Bmatrix = system.Bmatrix
for node in system.nodes:
if node.ideal_connected_with == '':
self.nodes.append(ResultsNode(topo_node=node))
......@@ -102,16 +102,18 @@ class Results():
def calculateI(self):
"""
To calculate the branch currents
Note: branch current flowing into start node coming from end node
"""
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_pu = - (self.nodes[fr].voltage_pu - self.nodes[to].voltage_pu) * self.Ymatrix[fr][to] + 1j*self.Bmatrix[fr][to] * self.nodes[fr].voltage_pu
branch.current = branch.current_pu * branch.topology_branch.base_current
def calculateIinj(self):
"""
calculate current injections at a node
Calculate current injections at a node
Note: node current flowing into the node
"""
for node in self.nodes:
to = complex(0, 0) # sum of the currents flowing to the node
......@@ -121,7 +123,7 @@ class Results():
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_pu = fr - to
node.current = node.current_pu * node.topology_node.base_current
def calculateSinj(self):
......@@ -152,7 +154,7 @@ class Results():
for node in self.nodes:
if branch_index == node.topology_node.index:
branch.power2_pu = -node.voltage_pu * (np.conj(branch.current_pu))
branch.power = branch.power2_pu * branch.topology_branch.base_apparent_power
branch.power2 = branch.power2_pu * branch.topology_branch.base_apparent_power
def get_node(self, index=None, uuid=None):
"""
......@@ -191,6 +193,22 @@ class Results():
voltages[node.topology_node.index] = node.voltage
return voltages
def get_branch_powers(self, pu=True):
"""
get branch powers
- if pu==True --> branch powers are expressed as per-unit
"""
#branch_powers = np.zeros(len(self.branches), dtype=np.complex_)
branch_powers = []
if pu == True:
for branch in self.branches:
branch_powers.append(branch.power_pu)
elif pu == False:
for branch in self.branches:
branch_powers.append(branch.power)
return branch_powers
def get_Iinj(self, pu=True):
"""
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment