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 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
from power_grid_model import 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.
idfield is a unique field across theinput_datadictionary. The keys ofinput-dataare predefined.
Please refer to Components for detailed explanation of all component types and their input/output attributes.
# node
node = initialize_array("input", "node", 3)
node["id"] = np.array([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]
# all
input_data = {
"node": node,
"line": line,
"sym_load": sym_load,
"source": source
}
We can print the input dataset by converting the numpy array to dataframe.
print(pd.DataFrame(input_data["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 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["node"]))
print("------line result------")
print(pd.DataFrame(output_data["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 component types in the 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 node.
output_data = model.calculate_power_flow(
symmetric=True,
error_tolerance=1e-8,
max_iterations=20,
calculation_method=CalculationMethod.newton_raphson,
output_component_types=["node"])
print("List of component types in result dataset")
print(list(output_data.keys()))
print("------node result------")
print(pd.DataFrame(output_data["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
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_toleranceandmax_iterationsas 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_toleranceandmax_iterationsare 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["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 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("update", "sym_load", 2)
update_sym_load["id"] = [4, 7] # same ID
update_sym_load["p_specified"] = [30e6, 15e6] # change active power
# leave reactive power the same, no need to specify
update_line = initialize_array("update", "line", 1)
update_line["id"] = [3] # change line ID 3
update_line["from_status"] = [0] # switch off at from side
# leave to-side swiching status the same, no need to specify
update_data = {
"sym_load": update_sym_load,
"line": update_line
}
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.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["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(
"update", "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["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["p_specified"] = [[30e6, 15e6]] * np.linspace(0, 1, 10).reshape(-1, 1)
time_series_mutation = {"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["line"]["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, ie. \(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["line"][0]))
id energized loading p_from q_from i_from \
0 3 1 0.212031 29511.239220 3.508895e+06 193.061627
1 5 1 0.064924 -1267.160636 -1.172098e+06 64.923606
2 8 1 0.156578 19805.345039 2.500724e+06 137.590869
s_from p_to q_to i_to s_to
0 3.509019e+06 1267.160636 -3.827902e+06 212.030847 3.827902e+06
1 1.172099e+06 3574.919119 8.320006e+05 46.000611 8.320083e+05
2 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["line"]["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("update", "line", (3, 3)) # 3 scenarios, 3 objects (lines)
# below the same broadcasting trick
line_profile["id"] = [[3, 5, 8]]
# fully specify the status of all lines, even it is the same as the base scenario
line_profile["from_status"] = [
[0, 1, 1],
[1, 0, 1],
[1, 1, 0]
]
line_profile["to_status"] = [
[0, 1, 1],
[1, 0, 1],
[1, 1, 0]
]
n_min_1_mutation_update_all = {"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["line"]["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("update", "line", (3, 1)) # 3 scenarios, 1 object mutation per scenario
# for each mutation, only one object is specified
line_profile["id"] = [[3], [5], [8]]
# specify only the changed status (switch off) of one line
line_profile["from_status"] = [[0], [0], [0]]
line_profile["to_status"] = [[0], [0], [0]]
n_min_1_mutation_update_specific = {"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["line"]["i_from"])
[[ 0. 1352.02947002 1962.69857764]
[1199.97577809 0. 573.32693369]
[1877.3560102 634.81512055 0. ]]
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 using exceptions in all cases.
Raise Exception
If there is an error in the ctypes wrapper, a relevant exception will be raised (e.g. KeyError). If there is an error in the C++ core, the C++ exception will be translated into a PowerGridError in Python. Below we show some examples for the errors in the construction, update, and calculation.
Construction
We try to construct a grid with a line connecting two nodes with different rated voltages. An error is raised in the construction.
from power_grid_model.errors import PowerGridError
# node
node_error = initialize_array("input", "node", 2)
node_error["id"] = [1, 2]
node_error["u_rated"] = [10.5e3, 150.0e3] # different rated voltages
# line
line_error = initialize_array("input", "line", 1)
line_error["id"] = [3]
line_error["from_node"] = [1]
line_error["to_node"] = [2]
line_error["from_status"] = [1]
line_error["to_status"] = [1]
line_error["r1"] = [0.25]
line_error["x1"] = [0.2]
line_error["c1"] = [10e-6]
line_error["tan1"] = [0.0]
line_error["i_n"] = [1000]
try:
PowerGridModel(
{
"node": node_error,
"line": line_error
}
)
except PowerGridError as e:
print(e)
Conflicting voltage for line 3
voltage at from node 1 is 10500.000000
voltage at to node 2 is 150000.000000
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.
line_update_error = initialize_array("update", "line", 1)
line_update_error["id"] = [12345] # non-existing
line_update_error["from_status"] = [1]
try:
model.update(update_data={"line": line_update_error})
except PowerGridError 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: 0.000000, error tolerance: 0.000000.
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["sym_load"]["id"][3] = 1000 # unknown id
time_series_mutation["sym_load"]["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', 'Iteration failed to converge after 20 iterations! Max deviation: 43691065954360128432752954469939668333328459052496602949289781148532517437499693334528.000000, error tolerance: 0.000000.\n']
# 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["node"]["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["node"]["u_pu"][e.succeeded_scenarios])
Node data with invalid results
[[0.99940117 0.99268579 0.99452137]
[0.99934769 0.98622639 0.98935286]
[0.99928838 0.97965401 0.98409554]
[0. 0. 0. ]
[0.99915138 0.96614948 0.97329879]
[0.99907317 0.95920586 0.96775071]
[0.9989881 0.95212621 0.96209647]
[0. 0. 0. ]
[0.99879613 0.93753005 0.95044796]
[0.9986885 0.92999747 0.94444167]]
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]]