test_nv_state_estimator.py 4.33 KB
Newer Older
Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
1
import os
2
import logging
Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
3
import numpy as np
4
5
6
import matplotlib.pyplot as plt

import cimpy
Markus Mirz's avatar
Markus Mirz committed
7
8
9
10
from pyvolt import network
from pyvolt import nv_state_estimator
from pyvolt import measurement
from pyvolt import results
11

12
logging.basicConfig(filename='CIGRE.log', level=logging.INFO, filemode='w')
13

14
xml_path = r"..\quickstart\sample_data\CIGRE-MV-NoTap"
Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
15
16
17
18
xml_files = [xml_path + r"\Rootnet_FULL_NE_06J16h_DI.xml",
                 xml_path + r"\Rootnet_FULL_NE_06J16h_EQ.xml",
                 xml_path + r"\Rootnet_FULL_NE_06J16h_SV.xml",
                 xml_path + r"\Rootnet_FULL_NE_06J16h_TP.xml"]
19

Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
20
21
22
23
xml_files_abs = []
for file in xml_files:
    xml_files_abs.append(os.path.abspath(file))
	
Jan Dinkelbach's avatar
Jan Dinkelbach committed
24
# read cim files and create new network.Systen object
Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
25
res, _ = cimpy.cim_import(xml_files_abs, "cgmes_v2_4_15")
26
system = network.System()
Jan Dinkelbach's avatar
Jan Dinkelbach committed
27
base_apparent_power = 25  # MW
Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
28
system.load_cim_data(res, base_apparent_power)
29

Jan Dinkelbach's avatar
Jan Dinkelbach committed
30
# read Input-Ergebnisdatei and store it in a results.Results object
31
loadflow_results_path = r"..\quickstart\sample_data\CIGRE-MV-NoTap"
Jan Dinkelbach's avatar
Jan Dinkelbach committed
32
loadflow_results_file = loadflow_results_path + r"\CIGRE-MV-NoTap.csv"
33
34
powerflow_results = results.Results(system)
powerflow_results.read_data_dpsim(loadflow_results_file)
35
36

# --- State Estimation with Ideal Measurements ---
Jan Dinkelbach's avatar
Jan Dinkelbach committed
37
""" Write here the percent uncertainties of the measurements"""
38
39
40
41
42
43
44
V_unc = 0
I_unc = 0
Sinj_unc = 0
S_unc = 0
Pmu_mag_unc = 0
Pmu_phase_unc = 0

45
46
# Create measurements object
"""use all node voltages as measures"""
Jan Dinkelbach's avatar
Jan Dinkelbach committed
47
measurements_set = measurement.MeasurementSet()
48
for node in powerflow_results.nodes:
Jan Dinkelbach's avatar
Jan Dinkelbach committed
49
50
    measurements_set.create_measurement(node.topology_node, measurement.ElemType.Node, measurement.MeasType.Vpmu_mag,
                                        np.absolute(node.voltage_pu), Pmu_mag_unc)
51
for node in powerflow_results.nodes:
Jan Dinkelbach's avatar
Jan Dinkelbach committed
52
53
    measurements_set.create_measurement(node.topology_node, measurement.ElemType.Node, measurement.MeasType.Vpmu_phase,
                                        np.angle(node.voltage_pu), Pmu_phase_unc)
54
measurements_set.meas_creation()
55
56

# Perform state estimation
Jan Dinkelbach's avatar
Jan Dinkelbach committed
57
state_estimation_results_ideal = nv_state_estimator.DsseCall(system, measurements_set)
58
59

# Show numerical comparison
Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
60
Vest_ideal = state_estimation_results_ideal.get_voltages(pu=False)
Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
61
Vtrue = powerflow_results.get_voltages(pu=False)
62
print(Vest_ideal - Vtrue)
63
64

# --- State Estimation with Non-Ideal Measurements ---
Jan Dinkelbach's avatar
Jan Dinkelbach committed
65
""" Write here the percent uncertainties of the measurements"""
66
67
68
Pmu_mag_unc = 1

# Create measurements data structures
69
"""use all node voltages as measures"""
Jan Dinkelbach's avatar
Jan Dinkelbach committed
70
measurements_set = measurement.MeasurementSet()
71
for node in powerflow_results.nodes:
Jan Dinkelbach's avatar
Jan Dinkelbach committed
72
73
    measurements_set.create_measurement(node.topology_node, measurement.ElemType.Node, measurement.MeasType.Vpmu_mag,
                                        np.absolute(node.voltage_pu), Pmu_mag_unc)
74
for node in powerflow_results.nodes:
Jan Dinkelbach's avatar
Jan Dinkelbach committed
75
76
    measurements_set.create_measurement(node.topology_node, measurement.ElemType.Node, measurement.MeasType.Vpmu_phase,
                                        np.angle(node.voltage_pu), Pmu_phase_unc)
77
measurements_set.meas_creation()
78
79

# Perform state estimation
Jan Dinkelbach's avatar
Jan Dinkelbach committed
80
state_estimation_results_real = nv_state_estimator.DsseCall(system, measurements_set)
81
82

# Show numerical comparison
Martin Alonso Moraga's avatar
Martin Alonso Moraga committed
83
Vest_real = state_estimation_results_real.get_voltages(pu=False)
84
print(Vest_real - Vtrue)
85
86

# Plot comparison
Jan Dinkelbach's avatar
Jan Dinkelbach committed
87
88
line_width = 6
fontsize = 26
89
90
91
92
93
94
95
96
97

plt.figure()
ax = plt.subplot()

nodes_uuid = [system.nodes[elem].uuid for elem in range(len(system.nodes))]

# Reorder and rescale results
idx_filter = np.argsort([int(uuid[1:]) for uuid in nodes_uuid])[1:]
nodes_uuid_filtered = [nodes_uuid[idx] for idx in idx_filter]
Jan Dinkelbach's avatar
Jan Dinkelbach committed
98
99
100
Vtrue_filtered = [abs(Vtrue[idx] / 1000) for idx in idx_filter]
Vest_ideal_filtered = [abs(Vest_ideal[idx] / 1000) for idx in idx_filter]
Vest_real_filtered = [abs(Vest_real[idx] / 1000) for idx in idx_filter]
101

Jan Dinkelbach's avatar
Jan Dinkelbach committed
102
plt.plot(Vest_ideal_filtered, linewidth=line_width, linestyle='-', label="state estimator (ideal measurements)")
103
plt.plot(Vtrue_filtered, linewidth=line_width, linestyle=':', label="DPsim load flow results")
Jan Dinkelbach's avatar
Jan Dinkelbach committed
104
plt.plot(Vest_real_filtered, linewidth=line_width, linestyle='-', label="state estimator (non-ideal measurements)")
105
106
107
108

plt.xticks(range(len(system.nodes)), fontsize=fontsize)
plt.yticks(fontsize=fontsize)
ax.set_xticklabels(nodes_uuid_filtered)
Jan Dinkelbach's avatar
Jan Dinkelbach committed
109
plt.ylabel("Node voltage [kV]", fontsize=fontsize)
Jan Dinkelbach's avatar
Jan Dinkelbach committed
110
plt.xlim([0, len(system.nodes) - 2])
Jan Dinkelbach's avatar
Jan Dinkelbach committed
111
plt.legend(fontsize=fontsize)
Jan Dinkelbach's avatar
Jan Dinkelbach committed
112
plt.show()