Power Flow Example

In this notebook we will walk through an example of power flow calculation of power-grid-model. The following points are covered

Construction of the model
Run power flow once, and its relevant function arguments
Update (change) of the model
Run power flow in batch calculations, and its relevant function arguments
Error handling

It serves as 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 (
    AttributeType,
    CalculationMethod,
    CalculationType,
    ComponentAttributeFilterOptions,
    ComponentType,
    DatasetType,
    LoadGenType,
    PowerGridModel,
    attribute_dtype,
    initialize_array,
)

Input Dataset

We create an input dataset by using the helper function initialize_array.

The units of all input are standard SI unit without any prefix, i.e. the unit of voltage is volt (V), not kV.
The attributes are of numpy array type.
AttributeType.id field is a unique field across the input_data dictionary. The keys of input-data are predefined.

Please refer to Components for detailed explanation of all component types and their input/output attributes.

# node
node = initialize_array(DatasetType.input, ComponentType.node, 3)
node[AttributeType.id] = np.array([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]

# all
input_data = {
    ComponentType.node: node,
    ComponentType.line: line,
    ComponentType.sym_load: sym_load,
    ComponentType.source: source,
}

It is also possible to specify a component input data in a columnar data format.

A columnar data format better integrates with most databases. In addition, it may bring memory and, in some cases, even computational performance improvements, because unused attribute columns can be omitted.

Make sure to provide the correct dtype to the numpy arrays, exposed for each dataset type, component and attribute via the helper function attribute_dtype function.

source_columns = {
    AttributeType.id: np.array([10], dtype=attribute_dtype(DatasetType.input, ComponentType.source, AttributeType.id)),
    AttributeType.node: np.array(
        [1], dtype=attribute_dtype(DatasetType.input, ComponentType.source, AttributeType.node)
    ),
    AttributeType.status: np.array(
        [1], dtype=attribute_dtype(DatasetType.input, ComponentType.source, AttributeType.status)
    ),
    AttributeType.u_ref: np.array(
        [1.0], dtype=attribute_dtype(DatasetType.input, ComponentType.source, AttributeType.u_ref)
    ),
    # We're not creating columns for u_ref_angle, sk, ... Instead, the default values are used. This saves us memory.
}

input_data = {
    ComponentType.node: node,
    ComponentType.line: line,
    ComponentType.sym_load: sym_load,
    ComponentType.source: source_columns,
}

We can print the input dataset by converting the numpy array to dataframe.

print(pd.DataFrame(input_data[ComponentType.sym_load]))

   id  node  status  type  p_specified  q_specified
0   4     2       1     0   20000000.0    5000000.0
1   7     6       1     0   10000000.0    2000000.0

(Assignment 1)

Validation (optional)

For efficiency reasons, most of the data is not explicitly validated in the power grid model. However, in most cases, a power flow calculation 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 the Error Handling section in this example notebook and 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.power_flow)

(Assignment 2)

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

(Assignment 3)

One-time Power Flow Calculation

You can call the method calculate_power_flow to do a one-time 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 power flow calculation with (they are also the default values for those arguments)

Symmetric calculation
Newton-Raphson method
Error tolerance: 1e-8
Maximum number of iteration: 20

output_data = model.calculate_power_flow(
    symmetric=True, error_tolerance=1e-8, max_iterations=20, calculation_method=CalculationMethod.newton_raphson
)

Result Dataset

We can also print a result dataset of node and line by converting the array to dataframe.

print("------node result------")
print(pd.DataFrame(output_data[ComponentType.node]))
print("------line result------")
print(pd.DataFrame(output_data[ComponentType.line]))

------node result------
   id  energized      u_pu             u   u_angle             p             q
0   1          1  0.998988  10489.375043 -0.003039  3.121451e+07  6.991358e+06
1   2          1  0.952126   9997.325181 -0.026031 -2.000000e+07 -5.000000e+06
2   6          1  0.962096  10102.012975 -0.021895 -1.000000e+07 -2.000000e+06
------line result------
   id  energized   loading        p_from        q_from      i_from  \
0   3          1  0.985666  1.736010e+07  4.072097e+06  981.460041   
1   5          1  0.205940 -3.365614e+06 -1.178649e+06  205.939917   
2   8          1  0.783206  1.385441e+07  2.919262e+06  779.311446   

         s_from          p_to          q_to        i_to          s_to  
0  1.783129e+07 -1.663439e+07 -3.821351e+06  985.666324  1.706768e+07  
1  3.566030e+06  3.396558e+06  8.861080e+05  200.617323  3.510241e+06  
2  1.415863e+07 -1.339656e+07 -2.886108e+06  783.206396  1.370392e+07  

Select components / attributes and format of result dataset

By default power-grid-model calculates the result attributes for all components. If you do not need all the component types in your results, for example if you are only interested in node, you can specify the list of the component types you want in the result dataset in the argument output_component_types.

The following code calculate the power flow and only generate the output attributes for ComponentType.node.

output_data = model.calculate_power_flow(
    symmetric=True,
    error_tolerance=1e-8,
    max_iterations=20,
    calculation_method=CalculationMethod.newton_raphson,
    output_component_types=[ComponentType.node],
)

print("List of component types in result dataset")
print(list(output_data.keys()))
print("------node result------")
print(pd.DataFrame(output_data[ComponentType.node]))

List of component types in result dataset
[node]
------node result------
   id  energized      u_pu             u   u_angle             p             q
0   1          1  0.998988  10489.375043 -0.003039  3.121451e+07  6.991358e+06
1   2          1  0.952126   9997.325181 -0.026031 -2.000000e+07 -5.000000e+06
2   6          1  0.962096  10102.012975 -0.021895 -1.000000e+07 -2.000000e+06

You can also request the output as a dictionary of attribute and their values columns for each component type.

output_data = model.calculate_power_flow(
    symmetric=True,
    error_tolerance=1e-8,
    max_iterations=20,
    calculation_method=CalculationMethod.newton_raphson,
    output_component_types={
        ComponentType.line: [AttributeType.id, AttributeType.p_from],  # line output columns id and p_from
    },
)
print("List of component types in result dataset")
print(list(output_data))
print("List of attribute types in line result")
print(list(output_data[ComponentType.line]))
print("------line result------")
print(output_data[ComponentType.line])

List of component types in result dataset
[line]
List of attribute types in line result
[id, p_from]
------line result------
{id: array([3, 5, 8], dtype=int32), p_from: array([17360100.20222374, -3365613.74450156, 13854413.52498137])}

You can also mix output types between components. In this example:

None requests row-based data
a list of attribute names requests columns for those attributes (as before)
ComponentAttributeFilterOptions.everything requests columns for all supported attributes of a component.

output_data = model.calculate_power_flow(
    symmetric=True,
    error_tolerance=1e-8,
    max_iterations=20,
    calculation_method=CalculationMethod.newton_raphson,
    output_component_types={
        # node output as row-based
        ComponentType.node: None,
        # line output columns id and p_from
        ComponentType.line: [AttributeType.id, AttributeType.p_from],
        # all sym_load attributes as columns
        ComponentType.sym_load: ComponentAttributeFilterOptions.everything,
    },
)

print("List of component types in result dataset")
print(list(output_data.keys()))
print("------node result------")
print(output_data[ComponentType.node].dtype.names)
print("------line result attributes------")
print(list(output_data[ComponentType.line].keys()))
print("------sym_load result attributes------")
print(list(output_data[ComponentType.sym_load].keys()))

List of component types in result dataset
[node, line, sym_load]
------node result------
(id, energized, u_pu, u, u_angle, p, q)
------line result attributes------
[id, p_from]
------sym_load result attributes------
[id, energized, p, q, i, s, pf]

Finally, it is possible to request columnar data for all attributes using the ComponentAttributeFilterOptions.everything as a shorthand.

output_data = model.calculate_power_flow(
    symmetric=True,
    error_tolerance=1e-8,
    max_iterations=20,
    calculation_method=CalculationMethod.newton_raphson,
    # all attributes for all component types as columns
    output_component_types=ComponentAttributeFilterOptions.everything,
)

print("List of component types in result dataset")
print(list(output_data.keys()))
print("------node result------")
print(list(output_data[ComponentType.node].keys()))
print("------line result attributes------")
print(list(output_data[ComponentType.line].keys()))
print("------sym_load result attributes------")
print(list(output_data[ComponentType.sym_load].keys()))

List of component types in result dataset
[node, line, sym_load, source]
------node result------
[id, energized, u_pu, u, u_angle, p, q]
------line result attributes------
[id, energized, loading, p_from, q_from, i_from, s_from, p_to, q_to, i_to, s_to]
------sym_load result attributes------
[id, energized, p, q, i, s, pf]

Calculation Method

There are four calculation methods available. They are listed in the CalculationMethod enumeration as follows:

Iterative methods: These Iterative methods converge to the true solution. You specify the error_tolerance and max_iterations as the iteration stopping condition. If the maximum number of iteration is reached and the error is still higher than the error tolerance, an exception will be thrown.
- CalculationMethod.newton_raphson: traditional Newton-Raphson method.
- CalculationMethod.iterative_current: a novel jacobi like method.
Linear methods: linearization approximation methods. In this case the error_tolerance and max_iterations are not relevant here. They are also garanteed to have a solution.
- CalculationMethod.linear: It treats all loads as constant impedances.
- CalculationMethod.linear_current: It is essentially a single iteration of Iterative current.

Following command calculates the power flow with linear method. Note the difference of result for node compared to the results from Newton-Raphson method.

output_data = model.calculate_power_flow(symmetric=True, calculation_method=CalculationMethod.linear)
print("------node result------")
print(pd.DataFrame(output_data[ComponentType.node]))

------node result------
   id  energized      u_pu             u   u_angle             p             q
0   1          1  0.999087  10490.414234 -0.002789  2.862536e+07  6.264884e+06
1   2          1  0.956325  10041.412991 -0.023873 -1.829115e+07 -4.572788e+06
2   6          1  0.965236  10134.973932 -0.020204 -9.316798e+06 -1.863360e+06

(Assignment 4)

Update Model

One can mutate the model by providing a dictionary of update data. This is similar to creating an input dataset. The only difference is to initialize the array with the DatasetType.update key.

Refer to Components for detailed explanation of all mutable attributes per component type.

NOTE: the initialized the array contains default NULL values for everything. If you do not need to update an attribute, you do not have to specify it. In the C++ core, the default NULL value will be treated as “not updating” for this attribute.

Create Update Dataset

The following code creates an update dataset which changes the two loads and switches Line #3 off at the from side.

update_sym_load = initialize_array(DatasetType.update, ComponentType.sym_load, 2)
update_sym_load[AttributeType.id] = [4, 7]  # same ID
update_sym_load[AttributeType.p_specified] = [30e6, 15e6]  # change active power
# leave reactive power the same, no need to specify

update_line = initialize_array(DatasetType.update, ComponentType.line, 1)
update_line[AttributeType.id] = [3]  # change line ID 3
update_line[AttributeType.from_status] = [0]  # switch off at from side
# leave to-side swiching status the same, no need to specify

update_data = {ComponentType.sym_load: update_sym_load, ComponentType.line: update_line}

If you are updating all components of the same type (uniformly) in the same order as in the input data, providing IDs is optional.

update_sym_load_no_id = initialize_array(DatasetType.update, ComponentType.sym_load, 2)
update_sym_load_no_id[AttributeType.p_specified] = [30e6, 15e6]  # change active power for both sym_loads
# leave reactive power the same, no need to specify

update_line_no_id = initialize_array(DatasetType.update, ComponentType.line, 3)
update_line_no_id[AttributeType.from_status] = [0, 1, 1]  # switch off at from side for line #3

update_data_no_id = {ComponentType.sym_load: update_sym_load_no_id, ComponentType.line: update_line_no_id}

Columnar data is also supported.

columnar_update_line = {
    AttributeType.id: np.array(
        [3], dtype=attribute_dtype(DatasetType.update, ComponentType.line, AttributeType.id)
    ),  # change line ID 3
    AttributeType.from_status: np.array(
        [0], dtype=attribute_dtype(DatasetType.update, ComponentType.line, AttributeType.from_status)
    ),  # switch off at from side
}
# leave to-side swiching status the same, no need to specify

update_data_col = {ComponentType.sym_load: update_sym_load, ComponentType.line: columnar_update_line}

Columnar data also supports optional IDs.

columnar_no_ID_update_line = {
    # Update IDs are not specified
    AttributeType.from_status: np.array(
        [0, 1, 1], dtype=attribute_dtype(DatasetType.update, ComponentType.line, AttributeType.from_status)
    ),  # The update for the whole column needs to be specified
}
# leave to-side swiching status the same, no need to specify

update_data_col_no_id = {ComponentType.sym_load: update_sym_load, ComponentType.line: columnar_no_ID_update_line}

Validation (optional)

When performing batch calculation/updating dataset, power-grid-model provides the functions to validate. See the Error Handling section in this example notebook, or Validation Examples for more detailed information on the validation functions.

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.power_flow)
assert_valid_batch_data(
    input_data=input_data, update_data=update_data_no_id, calculation_type=CalculationType.power_flow
)
assert_valid_batch_data(input_data=input_data, update_data=update_data_col, calculation_type=CalculationType.power_flow)
assert_valid_batch_data(
    input_data=input_data, update_data=update_data_col_no_id, calculation_type=CalculationType.power_flow
)

Call Update Method

Call the update function on a copy of the model. Note that you don’t need to make a copy in general, this is only done to keep the original model intact.

model_2 = model.copy()
model_2.update(update_data=update_data)

Re-calculate Power Flow

We calculate power flow again with all default arguments. And print the line results. See the difference of the updated data. Line #3 has zero power flow at the from-side and very small power flow at to-side (due to capacitance).

output_data = model_2.calculate_power_flow()
print("------line result------")
print(pd.DataFrame(output_data[ComponentType.line]))

------line result------
   id  energized   loading        p_from        q_from       i_from  \
0   3          1  0.014030  0.000000e+00  0.000000e+00     0.000000   
1   5          1  2.268159 -3.000004e+07 -4.812057e+06  2268.158655   
2   8          1  3.249928  5.676867e+07  1.571464e+07  3245.829906   

         s_from          p_to          q_to         i_to          s_to  
0  0.000000e+00  3.691991e+01 -1.879430e+05    14.030122  1.879430e+05  
1  3.038352e+07  3.385470e+07  7.678867e+06  2265.269969  3.471463e+07  
2  5.890358e+07 -4.885470e+07 -9.678867e+06  3249.927645  4.980424e+07  

Batch Calculation

We can use the same method calculate_power_flow to calculate a number of scenarios in one go. To do this, you need to supply an update_data argument. This argument contains a dictionary of 2D update arrays (one array per component type).

The model uses the current data as the base scenario. For each individual calculation, the model applies each mutation to the base scenario and calculates the power flow.

NOTE: after the batch calculation, the original model will be kept unchanged. Internally the program copies the original model to multiple batch models for the calculation.

Independent Batch Dataset

There are two ways to specify the mutations.

For each scenario only specify the objects that are changed in this scenario.
For each scenario specify all objects that are changed in one or more scenarios.

The latter is called independent batch dataset. Because all relevant objects are specified in each batch, different choices regarding performance optimization may be made in either case.

In general, the following is advised:

Use the non-independent batch dataset approach whenever few parameters change per scenario, but the batch samples many different components, e.g. during N-1 tests.
Use the independent batch dataset approach when a dense sampling of the parameter space is desired for relatively a few different components, e.g. during time series power flow calculation

See also performance guide for the latest recommendations.

Examples

Below we show three examples.

Time Series Profile

The following code creates a load profile with 10 timestamps for the two loads. The two loads are always present for all mutation scenarios.

load_profile = initialize_array(
    DatasetType.update, ComponentType.sym_load, (10, 2)
)  # note the shape of the array, 10 scenarios, 2 objects (loads)
# below is an assignment of shape (1, 2) array to shape (10, 2) array
# the numpy broadcasting rule ensures that the same object ids are repeated 10 times
# therefore the two objects are present for all the scenarios
load_profile[AttributeType.id] = [[4, 7]]
# this is a scale of load from 0% to 100%
# the array is an (10, 2) shape, each row is a scenario, each column is an object
load_profile[AttributeType.p_specified] = [[30e6, 15e6]] * np.linspace(0, 1, 10).reshape(-1, 1)

time_series_mutation = {ComponentType.sym_load: load_profile}

assert_valid_batch_data(
    input_data=input_data, update_data=time_series_mutation, calculation_type=CalculationType.power_flow
)

We can calculate the time series and print the current of the lines.

output_data = model.calculate_power_flow(update_data=time_series_mutation)
print(output_data[ComponentType.line][AttributeType.i_from])

[[ 193.06162675   64.92360593  137.59086941]
 [ 248.65360093   72.28087746  185.86691646]
 [ 368.12834615   90.73999329  285.46275837]
 [ 510.20016036  115.51167381  401.15982574]
 [ 662.04311447  143.73983633  523.57045909]
 [ 819.63118685  174.08971222  649.95104138]
 [ 981.46004118  205.93991655  779.31144601]
 [1146.90787963  238.98921914  911.25418138]
 [1315.7236725   273.08968816 1045.6266953 ]
 [1487.83778526  308.17354058 1182.39546494]]

Accessing batch data

It may be a bit unintuitive to read the output_data or update_data of a component directly as they are a dictionary of 3 dimension data, i.e., \(ids \times batches \times attributes\). Remember that the output_data or update_data are a dictionary of numpy structured arrays. Hence the component should be indexed first. The index that follows can be indexed with numpy structured arrays

To read the result of a single batch, lets say 1st batch,

print(pd.DataFrame(output_data[ComponentType.line][0]))

   id  energized   loading        p_from        q_from      i_from  \
 3          1  0.212031  29511.239220  3.508895e+06  193.061627   
 5          1  0.064924  -1267.160636 -1.172098e+06   64.923606   
 8          1  0.156578  19805.345039  2.500724e+06  137.590869   

         s_from         p_to          q_to        i_to          s_to  
3.509019e+06  1267.160636 -3.827902e+06  212.030847  3.827902e+06  
1.172099e+06  3574.919119  8.320006e+05   46.000611  8.320083e+05  
2.500802e+06 -3574.919119 -2.832001e+06  156.577603  2.832003e+06  

Or maybe we wish to find result of a single component, (e.g., 1st line) in all batches

print(output_data[ComponentType.line][AttributeType.i_from][:, 0])

[ 193.06162675  248.65360093  368.12834615  510.20016036  662.04311447
  819.63118685  981.46004118 1146.90787963 1315.7236725  1487.83778526]

(Assignment 5)

N-1 Scenario while specifying all objects in each scenario

The following code creates a N-1 scenario for all three lines. There are 3 scenarios, in each scenario, the from/to status of one line is switched off. However, all three lines are present in all mutation dataset. This means that for each scenario you will update two lines with the same switching status.

Specifying all objects can bring computational advantages.

line_profile = initialize_array(DatasetType.update, ComponentType.line, (3, 3))  # 3 scenarios, 3 objects (lines)
# below the same broadcasting trick
line_profile[AttributeType.id] = [[3, 5, 8]]
# fully specify the status of all lines, even it is the same as the base scenario
line_profile[AttributeType.from_status] = [[0, 1, 1], [1, 0, 1], [1, 1, 0]]
line_profile[AttributeType.to_status] = [[0, 1, 1], [1, 0, 1], [1, 1, 0]]

n_min_1_mutation_update_all = {ComponentType.line: line_profile}

assert_valid_batch_data(
    input_data=input_data, update_data=n_min_1_mutation_update_all, calculation_type=CalculationType.power_flow
)

We can calculate the N-1 and print the current of the lines. It is clear that per scenario one line is disabled.

output_data = model.calculate_power_flow(update_data=n_min_1_mutation_update_all)
print(output_data[ComponentType.line][AttributeType.i_from])

[[   0.         1352.02947002 1962.69857764]
 [1199.97577809    0.          573.32693369]
 [1877.3560102   634.81512055    0.        ]]

N-1 Scenario where only the changed objects are specified

The following code creates a N-1 scenario for all three lines. There are 3 scenarios, in each scenario, the from/to status of one line is switched off. In this dataset we only specify one line per mutation.

line_profile = initialize_array(
    DatasetType.update, ComponentType.line, (3, 1)
)  # 3 scenarios, 1 object mutation per scenario
# for each mutation, only one object is specified
line_profile[AttributeType.id] = [[3], [5], [8]]
# specify only the changed status (switch off) of one line
line_profile[AttributeType.from_status] = [[0], [0], [0]]
line_profile[AttributeType.to_status] = [[0], [0], [0]]

n_min_1_mutation_update_specific = {ComponentType.line: line_profile}

assert_valid_batch_data(
    input_data=input_data, update_data=n_min_1_mutation_update_specific, calculation_type=CalculationType.power_flow
)

We can calculate the N-1 and print the current of the lines. Note the results are the same as the independent N-1.

output_data = model.calculate_power_flow(update_data=n_min_1_mutation_update_specific)
print(output_data[ComponentType.line][AttributeType.i_from])

[[   0.         1352.02947002 1962.69857764]
 [1199.97577809    0.          573.32693369]
 [1877.3560102   634.81512055    0.        ]]

Cartesian product of batch datasets

It is possible to conduct a batch calculation of multiple datasets in form of a cartesian product of their scenarios. Assume certain batch datasets with N1, N2, N3, … scenarios. This would give us \(N1 * N2 * N3 * ...\) possible combinations via the cartesian product. The resultant output data is in flat form and it has dimension of N1 * N2 * N3 with first dataset being the highest dimension. This can be beneficial in reducing complexity of implementation of such batch calculation along with keeping the size of such resultant update_data to a minimum.

We can combine the time series mutation for all n-1 contingencies together in following way. Both the batch datasets are passed together in update_data in a list.

output_data = model.calculate_power_flow(update_data=[n_min_1_mutation_update_specific, time_series_mutation])
print("Output data has shape", output_data[ComponentType.line].shape)
line_output = output_data[ComponentType.line][AttributeType.energized]
print(line_output[0, :])
print(line_output[1, :])
print(line_output[10, :])
print(line_output[11, :])

Output data has shape (30, 3)
[0 1 1]
[0 1 1]
[1 0 1]
[1 0 1]

In addition, one can validate the cartesian product of batch datasets by validating its conforming datasets individually, as it is done in the validation example notebook. Conversely, the data validator in the validation module does not directly support a cartesian product of datasets, as this is not a new data structure, but just a list of datasets.

Caching Topology

If a batch scenario does not change the switching status of branches and sources, the model will re-use the pre-built internal graph/matrices for each calculation. Time-series load profile calculation is a typical use case. This can bring performance benefits.
If a batch scenario changes the switching status of branches and sources, the topology changes and no caching can be done for that scenario. N-1 check is a typical use case.
See the Performance Guide for tips and tricks regarding performance.

Parallel Computing

The batch calculation supports shared memory multi-threading parallel computing. The common internal states and variables are shared as much as possible to save memory usage and avoid copy. In the C++ implementation std::shared_ptr of const variable is used to share the resources between threads (which will not mutate the shared data).

You can set threading parameter to enable/disable parallel computing.

threading=-1, use sequential computing (default)
threading=0, use number of threads available from the machine hardware (recommended)
threading>0, set the number of threads you want to use

(Assignment 6)

Error Handling

The error handeling of power-grid-model is done with exceptions in all cases. We also provide an validation mechanism, which validates data structures and values offline. It is recommended to always validate your data before constructing a PowerGridModel instance. An alternative approach would be to validate only when an exception is raised, but be aware that not all data errors will raise exceptions: most of them will just yield invalid results without warning. Below we give examples of catching different error types in power grid model calculations. Refer to Validation Examples for more detailed information on the validation functions.

Construction

We try to construct a grid with a line connecting two nodes with different rated voltages. An ConflictVoltage error is caught in the construction.

from power_grid_model.errors import ConflictVoltage, PowerGridError

# node
node_error = initialize_array(DatasetType.input, ComponentType.node, 2)
node_error[AttributeType.id] = [1, 2]
node_error[AttributeType.u_rated] = [10.5e3, 150.0e3]  # different rated voltages
# line
line_error = initialize_array(DatasetType.input, ComponentType.line, 1)
line_error[AttributeType.id] = [3]
line_error[AttributeType.from_node] = [1]
line_error[AttributeType.to_node] = [2]
line_error[AttributeType.from_status] = [1]
line_error[AttributeType.to_status] = [1]
line_error[AttributeType.r1] = [0.25]
line_error[AttributeType.x1] = [0.2]
line_error[AttributeType.c1] = [10e-6]
line_error[AttributeType.tan1] = [0.0]
line_error[AttributeType.i_n] = [1000]

error_data = {ComponentType.node: node_error, ComponentType.line: line_error}

try:
    assert_valid_input_data(error_data, symmetric=True)
    model = PowerGridModel(error_data)
    output_data = model.calculate_state_estimation(symmetric=True)
except ConflictVoltage as e:
    print(e)

Conflicting voltage for line 3
 voltage at from node 1 is 10500
 voltage at to node 2 is 150000

Try validate_input_data() or validate_batch_data() to validate your data.

Update

We try to update the model with a non-existing line id.

from power_grid_model.errors import IDNotFound

line_update_error = initialize_array(DatasetType.update, ComponentType.line, 1)
line_update_error[AttributeType.id] = [12345]  # non-existing
line_update_error[AttributeType.from_status] = [1]


try:
    model.update(update_data={ComponentType.line: line_update_error})
except IDNotFound as e:
    print(e)

The id cannot be found: 12345

Try validate_input_data() or validate_batch_data() to validate your data.

Error in Calculation

Following command calculates the power flow with Newton-Raphson method that fails to converge, because the error tolerance is too low. It throws a PowerGridError

try:
    model.calculate_power_flow(
        symmetric=True, error_tolerance=1e-20, max_iterations=20, calculation_method=CalculationMethod.newton_raphson
    )
except PowerGridError as e:
    print(e)

Iteration failed to converge after 20 iterations! Max deviation: 3.54512293063893e-16, error tolerance: 1e-20.

Try validate_input_data() or validate_batch_data() to validate your data.

Error in Batch Calculation

Because parallel computing is allowed in batch calculation. The program will be terminated if an exception is thrown from the threads. Therefore, all the C++ exceptions are caught within the thread, and the error messages are re-thrown in the main thread and propagate as PowerGridBatchError in the Python side.

There could happen that only some of the scenarios fail in a batch calculation and others succeed. You can choose to specify the argument continue_on_batch_error (default False) if an error should be thrown in this case or let the program continue with partially valid results.

We modify the time-series load profile to add some errors. We add an unknown id for scenario 3, and a too large power in scenario 7.

time_series_mutation[ComponentType.sym_load][AttributeType.id][3] = 1000  # unknown id
time_series_mutation[ComponentType.sym_load][AttributeType.p_specified][7] = 1e100  # large power

# we run the batch calculation with default, it will raise error

from power_grid_model.errors import PowerGridBatchError

try:
    output_data = model.calculate_power_flow(
        update_data=time_series_mutation,
    )
except PowerGridBatchError as e:
    print(e)
    print(f"Failed scenarios: {e.failed_scenarios}")
    print(f"Succeeded scenarios: {e.succeeded_scenarios}")
    print(f"Error messages: {e.error_messages}")

There are errors in the batch calculation.
Try validate_input_data() or validate_batch_data() to validate your data.

Failed scenarios: [3 7]
Succeeded scenarios: [0 1 2 4 5 6 8 9]
Error messages: ['The id cannot be found: 1000\n', 'Sparse matrix error, possibly singular matrix!\nIf you get this error from state estimation, it might mean the system is not fully observable, i.e. not enough measurements.\nIt might also mean that you are running into a corner case where PGM cannot resolve yet.\nSee https://github.com/PowerGridModel/power-grid-model/issues/864.']

# we run the batch calculation with continue_on_batch_error=True,
#    it will return the results with partially valid data


output_data = model.calculate_power_flow(update_data=time_series_mutation, continue_on_batch_error=True)

# print node data for u_pu, note that the data is rubbish for scenario 3 and 7
print("Node data with invalid results")
print(output_data[ComponentType.node][AttributeType.u_pu])

# we can print the data with only succeeded scenarios
e = model.batch_error
print("Node data with only valid results")
print(output_data[ComponentType.node][AttributeType.u_pu][e.succeeded_scenarios])

Node data with invalid results
[[9.99401170e-001 9.92685785e-001 9.94521366e-001]
 [9.99347687e-001 9.86226389e-001 9.89352855e-001]
 [9.99288384e-001 9.79654011e-001 9.84095542e-001]
 [2.18565566e-312 4.31807050e-316 2.97079411e-313]
 [9.99151380e-001 9.66149483e-001 9.73298790e-001]
 [9.99073166e-001 9.59205860e-001 9.67750710e-001]
 [9.98988099e-001 9.52126208e-001 9.62096474e-001]
 [4.31810212e-316 4.31811082e-316 4.31810212e-316]
 [9.98796126e-001 9.37530046e-001 9.50447962e-001]
 [9.98688504e-001 9.29997471e-001 9.44441670e-001]]
Node data with only valid results
[[0.99940117 0.99268579 0.99452137]
 [0.99934769 0.98622639 0.98935286]
 [0.99928838 0.97965401 0.98409554]
 [0.99915138 0.96614948 0.97329879]
 [0.99907317 0.95920586 0.96775071]
 [0.9989881  0.95212621 0.96209647]
 [0.99879613 0.93753005 0.95044796]
 [0.9986885  0.92999747 0.94444167]]

Power flow calculations with automatic tap changing

For power flow calculations with automatic tap changing, please refer to the Transformer example notebook.