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 pandas as pd

from power_grid_model import (
    AngleMeasurementType,
    AttributeType,
    CalculationMethod,
    CalculationType,
    ComponentType,
    DatasetType,
    LoadGenType,
    MeasuredTerminalType,
    PowerGridModel,
    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 to Components for a detailed explanation of all component types and their input/output attributes.

# node
node = initialize_array(DatasetType.input, ComponentType.node, 3)
node[AttributeType.id] = [1, 2, 6]
node[AttributeType.u_rated] = [10.5e3, 10.5e3, 10.5e3]

# line
line = initialize_array(DatasetType.input, ComponentType.line, 3)
line[AttributeType.id] = [3, 5, 8]
line[AttributeType.from_node] = [1, 2, 1]
line[AttributeType.to_node] = [2, 6, 6]
line[AttributeType.from_status] = [1, 1, 1]
line[AttributeType.to_status] = [1, 1, 1]
line[AttributeType.r1] = [0.25, 0.25, 0.25]
line[AttributeType.x1] = [0.2, 0.2, 0.2]
line[AttributeType.c1] = [10e-6, 10e-6, 10e-6]
line[AttributeType.tan1] = [0.0, 0.0, 0.0]
line[AttributeType.i_n] = [1000, 1000, 1000]

# load
sym_load = initialize_array(DatasetType.input, ComponentType.sym_load, 2)
sym_load[AttributeType.id] = [4, 7]
sym_load[AttributeType.node] = [2, 6]
sym_load[AttributeType.status] = [1, 1]
sym_load[AttributeType.type] = [LoadGenType.const_power, LoadGenType.const_power]
sym_load[AttributeType.p_specified] = [20e6, 10e6]
sym_load[AttributeType.q_specified] = [5e6, 2e6]

# source
source = initialize_array(DatasetType.input, ComponentType.source, 1)
source[AttributeType.id] = [10]
source[AttributeType.node] = [1]
source[AttributeType.status] = [1]
source[AttributeType.u_ref] = [1.0]

# voltage sensor
sym_voltage_sensor = initialize_array(DatasetType.input, ComponentType.sym_voltage_sensor, 3)
sym_voltage_sensor[AttributeType.id] = [11, 12, 13]
sym_voltage_sensor[AttributeType.measured_object] = [1, 2, 6]
sym_voltage_sensor[AttributeType.u_sigma] = [1.0, 1.0, 1.0]
sym_voltage_sensor[AttributeType.u_measured] = [10489.37, 9997.32, 10102.01]
sym_voltage_sensor[AttributeType.u_angle_measured] = [0.0, 0.0, 0.0]

# power sensor
sym_power_sensor = initialize_array(DatasetType.input, ComponentType.sym_power_sensor, 5)
sym_power_sensor[AttributeType.id] = [14, 15, 16, 17, 18]
sym_power_sensor[AttributeType.measured_object] = [3, 5, 8, 4, 6]
sym_power_sensor[AttributeType.measured_terminal_type] = [
    MeasuredTerminalType.branch_to,
    MeasuredTerminalType.branch_from,
    MeasuredTerminalType.branch_to,
    MeasuredTerminalType.load,
    MeasuredTerminalType.node,
]
sym_power_sensor[AttributeType.power_sigma] = [1.0e3, 1.0e3, 1.0e3, 1.0e3, 1.0e3]
sym_power_sensor[AttributeType.p_measured] = [-1.66e07, -3.36e06, -1.33e07, 20e6, -10e6]
sym_power_sensor[AttributeType.q_measured] = [-3.82e06, -1.17e06, -2.88e06, 5e6, -2e6]

# current sensor
sym_current_sensor = initialize_array(DatasetType.input, ComponentType.sym_current_sensor, 3)
sym_current_sensor[AttributeType.id] = [19, 20, 21]
sym_current_sensor[AttributeType.measured_object] = [3, 5, 8]
sym_current_sensor[AttributeType.measured_terminal_type] = [
    MeasuredTerminalType.branch_from,
    MeasuredTerminalType.branch_to,
    MeasuredTerminalType.branch_from,
]
sym_current_sensor[AttributeType.angle_measurement_type] = [
    AngleMeasurementType.global_angle,
    AngleMeasurementType.local_angle,
    AngleMeasurementType.local_angle,
]
sym_current_sensor[AttributeType.i_sigma] = [100, 100, 100]
sym_current_sensor[AttributeType.i_angle_sigma] = [0.1, 0.1, 0.1]
sym_current_sensor[AttributeType.i_measured] = [978.814, 806.486, 776.753]
sym_current_sensor[AttributeType.i_angle_measured] = [-0.20864859135578062, 1.368, 1.349]

# all
input_data = {
    ComponentType.node: node,
    ComponentType.line: line,
    ComponentType.sym_load: sym_load,
    ComponentType.source: source,
    ComponentType.sym_voltage_sensor: sym_voltage_sensor,
    ComponentType.sym_power_sensor: sym_power_sensor,
    ComponentType.sym_current_sensor: sym_current_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[ComponentType.node]))

print("Line result")
display(pd.DataFrame(output_data[ComponentType.line]))

print("Sym_load result")
display(pd.DataFrame(output_data[ComponentType.sym_load]))

Node result

	id	energized	u_pu	u	u_angle	p	q
0	1	1	0.998883	10488.272639	0.013952	3.120033e+07	6.755258e+06
1	2	1	0.952278	9998.918864	-0.009381	-2.002141e+07	-4.812473e+06
2	6	1	0.962198	10103.079560	-0.005095	-9.968972e+06	-1.955178e+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.983966	1.736193e+07	3.929896e+06	979.903666	1.780115e+07	-1.663862e+07	-3.681085e+06	983.966281	1.704095e+07
1	5	1	0.205962	-3.382792e+06	-1.131388e+06	205.961755	3.566977e+06	3.413781e+06	8.387981e+05	200.886407	3.515321e+06
2	8	1	0.781260	1.383840e+07	2.825362e+06	777.480115	1.412388e+07	-1.338275e+07	-2.793976e+06	781.259502	1.367130e+07

Sym_load result

	id	energized	p	q	i	s	pf
0	4	1	2.002141e+07	4.812473e+06	1188.989174	2.059167e+07	0.972306
1	7	1	9.968972e+06	1.955178e+06	580.539808	1.015889e+07	0.981305

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.

In the following, state estimation is performed with the Newton-Raphson method. Note the difference in the results for the nodes compared to the results from the 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[ComponentType.node]))

Node result

	id	energized	u_pu	u	u_angle	p	q
0	1	1	0.999171	10491.294313	0.013903	3.114259e+07	6.983679e+06
1	2	1	0.952402	10000.219061	-0.009023	-1.997281e+07	-4.996913e+06
2	6	1	0.962376	10104.949737	-0.004879	-9.960845e+06	-2.000371e+06

It is also possible to specify the standard deviations of active and reactive power measurements 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(DatasetType.input, ComponentType.sym_power_sensor, 5)
sym_power_sensor2[AttributeType.id] = [14, 15, 16, 17, 18]
sym_power_sensor2[AttributeType.measured_object] = [3, 5, 8, 4, 6]
sym_power_sensor2[AttributeType.measured_terminal_type] = [
    MeasuredTerminalType.branch_to,
    MeasuredTerminalType.branch_from,
    MeasuredTerminalType.branch_to,
    MeasuredTerminalType.load,
    MeasuredTerminalType.node,
]
sym_power_sensor2[AttributeType.p_measured] = [-1.66e07, -3.36e06, -1.33e07, 20e6, -10e6]
sym_power_sensor2[AttributeType.q_measured] = [-3.82e06, -1.17e06, -2.88e06, 5e6, -2e6]
# define different standard deviations for p and q measurements
sym_power_sensor2[AttributeType.p_sigma] = [1.0e3, 1.0e4, 1.0e3, 1.0e4, 1.0e3]
sym_power_sensor2[AttributeType.q_sigma] = [1.0e4, 1.0e3, 1.0e4, 1.0e3, 1.0e4]

# all
input_data2 = {
    ComponentType.node: node,
    ComponentType.line: line,
    ComponentType.sym_load: sym_load,
    ComponentType.source: source,
    ComponentType.sym_voltage_sensor: sym_voltage_sensor,
    ComponentType.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[ComponentType.node]))

Node result

	id	energized	u_pu	u	u_angle	p	q
0	1	1	0.999144	10491.009558	0.013841	3.109941e+07	7.028384e+06
1	2	1	0.952418	10000.387584	-0.008985	-1.991963e+07	-5.002841e+06
2	6	1	0.962317	10104.325522	-0.004856	-9.973746e+06	-2.041416e+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: AttributeType.u and AttributeType.u_angle.

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

Where

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

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(DatasetType.update, ComponentType.sym_voltage_sensor, 1)
# for each mutation, only one object is specified
sym_voltage_sensor_update[AttributeType.id] = 13
sym_voltage_sensor_update[AttributeType.u_sigma] = np.inf  # disable this sensor

sym_power_sensor_update = initialize_array(DatasetType.update, ComponentType.sym_power_sensor, 1)
sym_power_sensor_update[AttributeType.id] = 18
sym_power_sensor_update[AttributeType.power_sigma] = np.inf  # disable this sensor

update_data = {
    ComponentType.sym_voltage_sensor: sym_voltage_sensor_update,
    ComponentType.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[ComponentType.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.998939	10488.854697	0.011360	3.113644e+07	6.783581e+06
1	2	1	0.952378	9999.970147	-0.011969	-2.004965e+07	-4.759283e+06
2	6	1	0.962303	10104.186383	-0.007496	-9.880385e+06	-2.039703e+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(DatasetType.update, ComponentType.sym_voltage_sensor, 3)
sym_voltage_sensor_update[AttributeType.id] = sym_voltage_sensor[AttributeType.id]
sym_voltage_sensor_update[AttributeType.u_sigma] = np.inf  # disable all sensors

sym_power_sensor_update = initialize_array(DatasetType.update, ComponentType.sym_power_sensor, 5)
sym_power_sensor_update[AttributeType.id] = sym_power_sensor[AttributeType.id]
sym_power_sensor_update[AttributeType.power_sigma] = np.inf  # disable all sensors

sym_current_sensor_update = initialize_array(DatasetType.update, ComponentType.sym_current_sensor, 3)
sym_current_sensor_update[AttributeType.id] = sym_current_sensor[AttributeType.id]
sym_current_sensor_update[AttributeType.i_sigma] = np.inf  # disable all sensors

update_data = {
    ComponentType.sym_voltage_sensor: sym_voltage_sensor_update,
    ComponentType.sym_power_sensor: sym_power_sensor_update,
    ComponentType.sym_current_sensor: sym_current_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.
No voltage sensor found!


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(
    DatasetType.update, ComponentType.sym_voltage_sensor, (4, 3)
)  # 4 scenarios, 3 objects per scenario
# for each mutation, only one object is specified
sym_voltage_sensor_update[AttributeType.id] = [[11, 12, 13]] * 4
sym_voltage_sensor_update[AttributeType.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[AttributeType.u_sigma][2, 2] = np.inf  # disable the third sensor of the third scenario

sym_power_sensor_update = initialize_array(DatasetType.update, ComponentType.sym_power_sensor, (4, 1))
sym_power_sensor_update[AttributeType.id] = [18]
sym_power_sensor_update[AttributeType.power_sigma] = [
    [1.0e3],
    [1.0e3],
    [1.0e3],
    [np.inf],  # disables this sensor on the last scenario
]

update_data = {
    ComponentType.sym_voltage_sensor: sym_voltage_sensor_update,
    ComponentType.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[ComponentType.node][AttributeType.u])

[[10465.44337926  9977.00924576 10080.55297677]
 [10480.47976245  9986.25423471 10093.14551907]
 [10415.61048312  9919.84715207 10025.79649187]
 [10451.69932534  9957.73806984 10066.52669313]]