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.

\[n_{\text{measurements}} >= n_{\text{unknowns}}\]

Where

\[n_{\text{unknowns}} = 2 n_{\text{nodes}} \]

And

\[n_{\text{measurements}} = n_{\text{nodes\_with\_voltage\_sensor\_without\_angle}} + 2 n_{\text{nodes\_with\_voltage\_sensor\_with\_angle}} + 2 n_{\text{branches\_with\_power\_sensor}} + 2 n_{\text{nodes\_without\_any\_appliances\_connected}} + 2 n_{\text{nodes\_with\_all\_connected\_appliances\_measured\_by\_power\_sensor}}\]

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)

Not enough measurements available for state estimation.

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]]