State Estimation Example
In this notebook we will walk through an example of state estimation calculation of power-grid-model
.
The following points are covered
Construction of the model
Run state estimation once, and its relevant function arguments
It serves an example of how to use the Python API. For detailed API documentation, refer to Python API reference and Native Data Interface.
Example Network
We use a simple network with 3 nodes, 1 source, 3 lines, and 2 loads. As shown below:
-----------------------line_8---------------
| |
node_1 ---line_3--- node_2 ----line_5---- node_6
| | |
source_10 sym_load_4 sym_load_7
The 3 nodes are connected in a triangular way by 3 lines.
# some basic imports
import numpy as np
import warnings
with warnings.catch_warnings(action="ignore", category=DeprecationWarning):
# suppress warning about pyarrow as future required dependency
import pandas as pd
from power_grid_model import LoadGenType
from power_grid_model import PowerGridModel, CalculationMethod, CalculationType, MeasuredTerminalType
from power_grid_model import initialize_array
Input Dataset
We create input dataset by using the helper function initialize_array
.
Note the units of all input are standard SI unit without any prefix,
i.e. the unit of voltage is volt (V), not kV.
Please refer Components for detailed explanation of all component types and their input/output attributes.
# node
node = initialize_array("input", "node", 3)
node["id"] = [1, 2, 6]
node["u_rated"] = [10.5e3, 10.5e3, 10.5e3]
# line
line = initialize_array("input", "line", 3)
line["id"] = [3, 5, 8]
line["from_node"] = [1, 2, 1]
line["to_node"] = [2, 6, 6]
line["from_status"] = [1, 1, 1]
line["to_status"] = [1, 1, 1]
line["r1"] = [0.25, 0.25, 0.25]
line["x1"] = [0.2, 0.2, 0.2]
line["c1"] = [10e-6, 10e-6, 10e-6]
line["tan1"] = [0.0, 0.0, 0.0]
line["i_n"] = [1000, 1000, 1000]
# load
sym_load = initialize_array("input", "sym_load", 2)
sym_load["id"] = [4, 7]
sym_load["node"] = [2, 6]
sym_load["status"] = [1, 1]
sym_load["type"] = [LoadGenType.const_power, LoadGenType.const_power]
sym_load["p_specified"] = [20e6, 10e6]
sym_load["q_specified"] = [5e6, 2e6]
# source
source = initialize_array("input", "source", 1)
source["id"] = [10]
source["node"] = [1]
source["status"] = [1]
source["u_ref"] = [1.0]
# voltage sensor
sym_voltage_sensor = initialize_array("input", "sym_voltage_sensor", 3)
sym_voltage_sensor["id"] = [11, 12, 13]
sym_voltage_sensor["measured_object"] = [1, 2, 6]
sym_voltage_sensor["u_sigma"] = [1.0, 1.0, 1.0]
sym_voltage_sensor["u_measured"] = [10489.37, 9997.32, 10102.01]
# power sensor
sym_power_sensor = initialize_array("input", "sym_power_sensor", 8)
sym_power_sensor["id"] = [14, 15, 16, 17, 18, 19, 20, 21]
sym_power_sensor["measured_object"] = [3, 3, 5, 5, 8, 8, 4, 6]
sym_power_sensor["measured_terminal_type"] = [
MeasuredTerminalType.branch_from,
MeasuredTerminalType.branch_to,
MeasuredTerminalType.branch_from,
MeasuredTerminalType.branch_to,
MeasuredTerminalType.branch_from,
MeasuredTerminalType.branch_to,
MeasuredTerminalType.load,
MeasuredTerminalType.node,
]
sym_power_sensor["power_sigma"] = [1.0e3, 1.0e3, 1.0e3, 1.0e3, 1.0e3, 1.0e3, 1.0e3, 1.0e3]
sym_power_sensor["p_measured"] = [1.73e07, -1.66e07, -3.36e06, 3.39e06, 1.38e07, -1.33e07, 20e6, -10e6]
sym_power_sensor["q_measured"] = [4.07e06, -3.82e06, -1.17e06, 8.86e05, 2.91e06, -2.88e06, 5e6, -2e6]
# all
input_data = {
"node": node,
"line": line,
"sym_load": sym_load,
"source": source,
"sym_voltage_sensor": sym_voltage_sensor,
"sym_power_sensor": sym_power_sensor,
}
Validation (optional)
For efficiency reasons, most of the data is not explicitly validated in the power grid model. However, in most cases, a state estimation will fail/crash if the data is invalid. Often with a low level error message that is hard to relate to the original objects. Therfore, it is recommended to always validate your data before constructing a PowerGridModel instance.
The simplest and most effective way to validate your data is by using assert_valid_input_data()
which will throw an error if it encounters any invalid data. See Validation Examples for more detailed information on the validation functions.
from power_grid_model.validation import assert_valid_input_data
assert_valid_input_data(input_data=input_data, calculation_type=CalculationType.state_estimation)
Construction
The construction of the model is just calling the constructor of PowerGridModel
. The default system_frequency
is 50.0 Hz.
model = PowerGridModel(input_data)
# model = PowerGridModel(input_data, system_frequency=60.0) # if you have another system frequency
One-time State Estimation Calculation
You can call the method calculate_state_estimation
to do a one-time state estimation calculation based on the current network data in the model. In this case you should not specify the argument update_data
as it is used for batch calculation.
The following command executes a one-time state estimation calculation with (they are also the default values for those arguments)
Symmetric calculation
Iterative linear method
Error tolerance: 1e-8
Maximum number of iteration: 20.
output_data = model.calculate_state_estimation(
symmetric=True, error_tolerance=1e-8, max_iterations=20, calculation_method=CalculationMethod.iterative_linear
)
Result Dataset
We can also print one result dataset of node and line by converting the array to dataframe.
print("Node result")
display(pd.DataFrame(output_data["node"]))
print("Line result")
display(pd.DataFrame(output_data["line"]))
print("Sym_load result")
display(pd.DataFrame(output_data["sym_load"]))
Node result
id | energized | u_pu | u | u_angle | p | q | |
---|---|---|---|---|---|---|---|
0 | 1 | 1 | 0.998922 | 10488.678869 | 0.000000 | 3.112866e+07 | 6.984619e+06 |
1 | 2 | 1 | 0.952157 | 9997.648144 | -0.022921 | -1.996348e+07 | -4.998968e+06 |
2 | 6 | 1 | 0.962132 | 10102.384509 | -0.018780 | -9.956656e+06 | -1.999083e+06 |
Line result
id | energized | loading | p_from | q_from | i_from | s_from | p_to | q_to | i_to | s_to | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 3 | 1 | 0.983446 | 1.731816e+07 | 4.068528e+06 | 979.232855 | 1.778965e+07 | -1.659573e+07 | -3.820392e+06 | 983.446244 | 1.702978e+07 |
1 | 5 | 1 | 0.206048 | -3.367751e+06 | -1.178575e+06 | 206.048359 | 3.568023e+06 | 3.398729e+06 | 8.860393e+05 | 200.729011 | 3.512325e+06 |
2 | 8 | 1 | 0.780866 | 1.381050e+07 | 2.916091e+06 | 776.961713 | 1.411500e+07 | -1.335538e+07 | -2.885123e+06 | 780.865585 | 1.366346e+07 |
Sym_load result
id | energized | p | q | i | s | pf | |
---|---|---|---|---|---|---|---|
0 | 4 | 1 | 1.996348e+07 | 4.998968e+06 | 1188.457478 | 2.057985e+07 | 0.970050 |
1 | 7 | 1 | 9.956656e+06 | 1.999083e+06 | 580.377789 | 1.015536e+07 | 0.980434 |
Calculation Method
There are two calculation methods available. They are listed in the CalculationMethod
enumeration as follows:
Iterative linear method:
CalculationMethod.iterative_linear
: All measurements are linearized before running the method.
Newton Raphson method:
CalculationMethod.newton_raphson
: traditional Newton-Raphson method.
Following command does state estimation with Newton Raphson method. Note the difference of result for node compared to the results from iterative linear method.
output_data_NR = model.calculate_state_estimation(symmetric=True, calculation_method=CalculationMethod.newton_raphson)
print("Node result")
display(pd.DataFrame(output_data_NR["node"]))
Node result
id | energized | u_pu | u | u_angle | p | q | |
---|---|---|---|---|---|---|---|
0 | 1 | 1 | 0.999187 | 10491.465789 | 0.000000 | 3.112736e+07 | 6.981011e+06 |
1 | 2 | 1 | 0.952440 | 10000.623189 | -0.022911 | -1.996275e+07 | -4.997141e+06 |
2 | 6 | 1 | 0.962411 | 10105.311164 | -0.018773 | -9.956895e+06 | -1.998502e+06 |
It is also possible to specify the standard deviations of active and reactive power mesurements seperately. This is especially powerful when using the Newton Raphson method for state estimation. Refer to the Documentation for a detailed explanation. Below is an example of state estimation with seperate standard deviations for p and q measurements.
# power sensor
sym_power_sensor2 = initialize_array("input", "sym_power_sensor", 8)
sym_power_sensor2["id"] = [14, 15, 16, 17, 18, 19, 20, 21]
sym_power_sensor2["measured_object"] = [3, 3, 5, 5, 8, 8, 4, 6]
sym_power_sensor2["measured_terminal_type"] = [
MeasuredTerminalType.branch_from,
MeasuredTerminalType.branch_to,
MeasuredTerminalType.branch_from,
MeasuredTerminalType.branch_to,
MeasuredTerminalType.branch_from,
MeasuredTerminalType.branch_to,
MeasuredTerminalType.load,
MeasuredTerminalType.node,
]
sym_power_sensor2["p_measured"] = [1.73e07, -1.66e07, -3.36e06, 3.39e06, 1.38e07, -1.33e07, 20e6, -10e6]
sym_power_sensor2["q_measured"] = [4.07e06, -3.82e06, -1.17e06, 8.86e05, 2.91e06, -2.88e06, 5e6, -2e6]
sym_power_sensor2["p_sigma"] = [1.0e3, 1.0e3, 1.0e4, 1.0e3, 1.0e4, 1.0e3, 1.0e4, 1.0e3] # define different standard deviations for p and q measurements
sym_power_sensor2["q_sigma"] = [1.0e4, 1.0e4, 1.0e3, 1.0e4, 1.0e3, 1.0e4, 1.0e3, 1.0e4]
# all
input_data2 = {
"node": node,
"line": line,
"sym_load": sym_load,
"source": source,
"sym_voltage_sensor": sym_voltage_sensor,
"sym_power_sensor": sym_power_sensor2,
}
model2 = PowerGridModel(input_data2)
output_data_NR2 = model2.calculate_state_estimation(symmetric=True, calculation_method=CalculationMethod.newton_raphson)
print("Node result")
display(pd.DataFrame(output_data_NR2["node"]))
Node result
id | energized | u_pu | u | u_angle | p | q | |
---|---|---|---|---|---|---|---|
0 | 1 | 1 | 0.998937 | 10488.835308 | 0.000000 | 3.110692e+07 | 6.979313e+06 |
1 | 2 | 1 | 0.952221 | 9998.324909 | -0.022891 | -1.993051e+07 | -4.995857e+06 |
2 | 6 | 1 | 0.962157 | 10102.645852 | -0.018779 | -9.969948e+06 | -1.998607e+06 |
Observability
In order to perform a state estimation the number of measurements should be larger than or equal to the number of unknowns. For each node there are two unknowns: u
and u_angle
.
Where
And
In this formula, only measurements with non-negligible contributions should be considered. That means that measurements with infinite variances have vanishing contribution and are to be omitted, regardless of whether those divergences are due to infinite sigmas or unreasonably large ones. As a result, missing data can be representing by providing infinity as a standard deviation, and the underdetermined edge case is handled correctly.
Some missing measurements
Few missing measurements may keep the data observable.
sym_voltage_sensor_update = initialize_array("update", "sym_voltage_sensor", 1)
# for each mutation, only one object is specified
sym_voltage_sensor_update["id"] = 13
sym_voltage_sensor_update["u_sigma"] = np.inf # disable this sensor
sym_power_sensor_update = initialize_array("update", "sym_power_sensor", 1)
sym_power_sensor_update["id"] = 21
sym_power_sensor_update["power_sigma"] = np.inf # disable this sensor
update_data = {"sym_voltage_sensor": sym_voltage_sensor_update, "sym_power_sensor": sym_power_sensor_update}
# Permanent update
model.update(update_data=update_data)
output_data = model.calculate_state_estimation()
display(pd.DataFrame(output_data["node"]))
# fully reinstigate the original model to prevent this cell from influencing re-runs. A model update would be sufficient as well.
model = PowerGridModel(input_data)
id | energized | u_pu | u | u_angle | p | q | |
---|---|---|---|---|---|---|---|
0 | 1 | 1 | 0.998927 | 10488.737113 | 3.391770e-21 | 3.110041e+07 | 6.980842e+06 |
1 | 2 | 1 | 0.952178 | 9997.866792 | -2.291293e-02 | -1.997486e+07 | -4.998422e+06 |
2 | 6 | 1 | 0.962193 | 10103.023556 | -1.874179e-02 | -9.918677e+06 | -1.997233e+06 |
Unobservable case
If too many measurements are invalid, the result may become unobservable.
from power_grid_model.errors import PowerGridError
sym_voltage_sensor_update = initialize_array("update", "sym_voltage_sensor", 3)
sym_voltage_sensor_update["id"] = sym_voltage_sensor["id"]
sym_voltage_sensor_update["u_sigma"] = np.inf # disable all sensors
sym_power_sensor_update = initialize_array("update", "sym_power_sensor", 8)
sym_power_sensor_update["id"] = sym_power_sensor["id"]
sym_power_sensor_update["power_sigma"] = np.inf # disable all sensors
update_data = {"sym_voltage_sensor": sym_voltage_sensor_update, "sym_power_sensor": sym_power_sensor_update}
try:
# Permanent update
model.update(update_data=update_data)
output_data = model.calculate_state_estimation()
print("Success")
except PowerGridError as e:
print(e)
finally:
# fully reinstigate the original model to prevent this cell from influencing re-runs. A model update would be sufficient as well.
model = PowerGridModel(input_data)
Sparse matrix error, possibly singular matrix!
If you get this error from state estimation, it usually means the system is not fully observable, i.e. not enough measurements.
Try validate_input_data() or validate_batch_data() to validate your data.
Batch calculation
For state estimation a batch calculation can be performed in a similar manner as for powerflow calculations. Below is a short example provided for reference. For more a more detailed example about performing a batch calculation, please refer to the Power Flow Example or visit the Documentation.
Create Update Dataset
The following code creates an update dataset which runs state estimation for a batch with three different measurements for the symmetric voltage sensors
sym_voltage_sensor_update = initialize_array(
"update", "sym_voltage_sensor", (4, 3)
) # 4 scenarios, 3 objects per scenario
# for each mutation, only one object is specified
sym_voltage_sensor_update["id"] = [[11, 12, 13]] * 4
sym_voltage_sensor_update["u_measured"] = [
[10402.31, 10023.68, 10095.40],
[10583.20, 9869.85, 10105.28],
[10482.87, 9850.85, 10140.77],
[10482.87, 9850.85, 10140.77],
]
sym_voltage_sensor_update["u_sigma"][2, 2] = np.inf # disable the third sensor of the third scenario
sym_power_sensor_update = initialize_array("update", "sym_power_sensor", (4, 1))
sym_power_sensor_update["id"] = [21]
sym_power_sensor_update["power_sigma"] = [
[1.0e3],
[1.0e3],
[1.0e3],
[np.inf], # disables this sensor on the last scenario
]
update_data = {"sym_voltage_sensor": sym_voltage_sensor_update, "sym_power_sensor": sym_power_sensor_update}
Validation (optional)
from power_grid_model.validation import assert_valid_batch_data
assert_valid_batch_data(
input_data=input_data, update_data=update_data, calculation_type=CalculationType.state_estimation
)
Call Update Method
batch_output_data = model.calculate_state_estimation(
update_data=update_data,
symmetric=True,
error_tolerance=1e-8,
max_iterations=20,
calculation_method=CalculationMethod.iterative_linear,
)
print(batch_output_data["node"]["u"])
[[10466.15445067 9974.61085414 10078.90372425]
[10480.594761 9987.299152 10094.19520541]
[10416.3953498 9920.53736882 10027.28386252]
[10452.49360638 9957.85948757 10066.06548575]]