We can define N-1 cases to be analyzed as contingencies. This means that indices of net.line, net.trafo, net.trafo3w can be defined as contingencies, which are switched off one at a time. The power system is analyzed with power flow calculations, and the min/max values among all the N-1 cases are obtained for relevant variables.
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import pandapower as pp
import pandapower.networks
import pandapower.contingency
import pandapower.control
import pandapower.timeseries
import pandapower.plotting
A dictionary that defines all the relevant N-1 cases to be analyzed is a required input parameter. Further parameters are:
The results of the contingency analysis are returned as a dictionary, and also written in the results tables in net:
net = pp.networks.case9()
nminus1_cases = {"line": {"index": net.line.index.values}}
res = pp.contingency.run_contingency(net, nminus1_cases)
The results can be evaluated with the fuction check_elements_within_limits.
It uses the required input element_limits, which is a dictionary with limit values for all the relevant limits. Elements with at least one limit variable set as not nan are considered relevant, and elements with all limit variables set to nan are ignored.
The function returns True if all results are within specified limits, and False in case of any violations:
element_limits = pp.contingency.get_element_limits(net)
pp.contingency.check_elements_within_limits(element_limits, res, True)
False
The contingency analysis determined violations of limits. The function report_contingency_results can be used to print logging statements, informing the user of the details of the limit violations:
pp.contingency.report_contingency_results(element_limits, res)
hp.pandapower.contingency.contingency - INFO: bus: min_vm_pu violation at index [8] ([0.794])
The function run_contingency can be passed to run_timeseries.
The parameters for run_contingency can be passed as kwargs to the run_timeseries function as the parameter run_control_fct. Note: avoid passing it as the run argument, because it will cause multiple contingency calculations per time step if there are any controllers in the grid model. Here we pass nminus1_cases as a keyword argument, which is passed to run_contingency during time-series simulation.
This way, the contingency analysis is executed at every timne step of the time-series simulation. The results of the contingency analysis for every time step can be added to the output of the output_writer and analyzed after the time-series calculation.
First, we configure the time-series simulation:
load_profiles = pd.DataFrame(net.load.p_mw.values * (np.random.random((24, len(net.load))) * 0.4 + 0.8),
index=np.arange(24), columns=net.load.index.values)
dsl = pp.timeseries.DFData(load_profiles)
pp.control.ConstControl(net, element="load", variable="p_mw", element_index=net.load.index.values,
profile_name=net.load.index.values, data_source=dsl)
gen_profiles = pd.DataFrame(net.gen.p_mw.values * (np.random.random((24, len(net.gen))) * 0.4 + 0.8),
index=np.arange(24), columns=net.gen.index.values)
dsg = pp.timeseries.DFData(gen_profiles)
pp.control.ConstControl(net, element="gen", variable="p_mw", element_index=net.gen.index.values,
profile_name=net.gen.index.values, data_source=dsg)
ow = pp.timeseries.OutputWriter(net)
ow.log_variable("res_line", "max_loading_percent")
ow.log_variable("res_line", "min_loading_percent")
ow.log_variable("res_bus", "max_vm_pu")
ow.log_variable("res_bus", "min_vm_pu")
The time-series simulation is executed:
pp.timeseries.run_timeseries(net, time_steps=range(24),
run_control_fct=pp.contingency.run_contingency, nminus1_cases=nminus1_cases)
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 24/24 [00:02<00:00, 10.70it/s]
We can aggregate the result values from the output writer and show the max. and min. values for the N-0 and N-1 situations for every time step:
fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(10,4))
ow.output["res_line.max_loading_percent"].max(axis=1).plot(label="max. Loading, N-1 (%)", ax=ax1, linestyle="--")
ow.output["res_line.loading_percent"].max(axis=1).plot(label="max. Loading, N-0 (%)", ax=ax1)
ow.output["res_line.min_loading_percent"].min(axis=1).plot(label="min. Loading, N-1 (%)", ax=ax1, linestyle="--")
ow.output["res_line.loading_percent"].min(axis=1).plot(label="min. Loading, N-0 (%)", ax=ax1)
ax1.legend(ncol=2, loc="upper center", bbox_to_anchor=[0.5, -0.15])
ax1.set_ylabel("Line loading (%)")
ax1.set_xlabel("Time step")
ow.output["res_bus.max_vm_pu"].max(axis=1).plot(label="max. $V_m$, N-1 (p.u.)", ax=ax2, linestyle="--")
ow.output["res_bus.vm_pu"].max(axis=1).plot(label="max. $V_m$, N-0 (p.u.)", ax=ax2)
ow.output["res_bus.min_vm_pu"].min(axis=1).plot(label="min. $V_m$, N-1 (p.u.)", ax=ax2, linestyle="--")
ow.output["res_bus.vm_pu"].min(axis=1).plot(label="min. $V_m$, N-0 (p.u.)", ax=ax2)
ax2.legend(ncol=2, loc="upper center", bbox_to_anchor=[0.5, -0.15])
ax2.set_ylabel("Bus voltage magnitude (p.u.)")
ax2.set_xlabel("Time step")
plt.tight_layout()
A more advanced plotting of results using boxplot provides more insight into every time step:
def plot_results_boxplots():
fig, ((ax11, ax12), (ax21, ax22)) = plt.subplots(ncols=2, nrows=2, sharey="row", sharex="col", figsize=(10,6))
ow.output["res_line.max_loading_percent"].T.boxplot(ax=ax11, whis=[0,100], grid=False)
ax11.set_title("max. Loading of lines per time step, N-1 (%)")
ax11.set_ylabel("Line loading (%)")
ow.output["res_line.min_loading_percent"].T.boxplot(ax=ax12, whis=[0,100], grid=False)
ax12.set_title("min. Loading of lines per time step, N-1 (%)")
ow.output["res_bus.max_vm_pu"].T.boxplot(ax=ax21, whis=[0,100], grid=False)
ax21.set_title("max. $V_m$ of buses per time step, N-1 (p.u.)")
ax21.set_ylabel("Bus voltage magnitude (p.u.)")
ax21.set_xlabel("Time step")
ow.output["res_bus.min_vm_pu"].T.boxplot(ax=ax22, whis=[0,100], grid=False)
ax22.set_title("min. $V_m$ of buses per time step, N-1 (p.u.)")
ax22.set_xlabel("Time step")
plt.tight_layout()
The library lightsim2grid implements its own contingency analysis function in C++, which provides at least a magnitude of performance improvement in comparison to the pandapower function. The caveats are that not all possible features (e.g. trafo3w) of pandapower grid modeling are compatible with this particular implementation, and that buses must have a continuous index starting with 0. If the pandapower grid model has incompatible features, a warning is shown.
Furthermore, there can be a difference in results between pandapower and lightsim2grid calculations due to some implementation differences. The results of the contingency analysis are related to the "from_bus" of lines and the "hv_bus" of transformers in lightsim2grid, but in pandapower the maximum value between "from_bus" and "to_bus" ("hv_bus" and "lv_bus") is used. In addition, if a contingency leads to an isolated island in the grid, the results for the contingency are set to nan in lightsim2grid. However, pandapower selects a gen element authomatically in such an isolated island as a new slack, and has values as results for such a contingency.
If the pandapower grid model has an ext_grid element as slack (as is most commonly the case), lightsim2grid will show warnings that it will attempt to use it as a slack generator. In order to avoid such warnings, we advise converting the ext_grid element to a gen element with True for slack parameter and set the slack_weight parameter to 1. This can be accomplished using the convenience function pandapower.toolbox.replace_ext_grid_by_gen()
.
It is important to know about these differences and consider them in configuring the calculation and in the analysis of the results. In the following example, we demonstrate the use of this function.
gen_index = pp.replace_ext_grid_by_gen(net, slack=True, cols_to_keep=["slack_weight"])
pp.timeseries.run_timeseries(net, time_steps=range(24),
run_control_fct=pp.contingency.run_contingency_ls2g, nminus1_cases=nminus1_cases)
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 24/24 [00:00<00:00, 62.27it/s]
The performance improvement is evident from comparing the number of iterations per second between the time-series simulations. We can observe an improvement of about 6 times here. However, a much higher improvement can be expected with larger grids.
The difference in results is obvious from looking at the plots below. This is not due to any errors but due to the aforementioned differences in the implementation. We strongly advise knowing exactly what those differences are when using this function.
plot_results_boxplots()
In addition to the impact of contingencies on the grid, it is also interesting to include the focus on the causes of the contingencies. The columns causes_overloading, cause_element, and cause_index describe which element ("line" or "trafo") caused the highest loading of the line. In the following example we color the lines according to the metric of how often they cause the highest overloading of other lines. This is useful to identify the most critical lines. The function run_contingency_ls2g additionally outputs the overall grid congestion caused by the failure of each elements i.e. for every N-1 case, which is useful to evaluate how critical any element is to the overall grid.
pp.plotting.set_line_geodata_from_bus_geodata(net)
net.load.p_mw = [125, 150, 175]
pp.contingency.run_contingency_ls2g(net, nminus1_cases)
pp.plotting.convert_geodata_to_gis(net)
d = pd.DataFrame(net.res_line.groupby("cause_index")["cause_index"].count())
net.line["cause_number"] = 0
net.line.loc[d.index, "cause_number"] = d.cause_index.values
Three lines are responsible for the highest loadings of other lines, and the line 5 has that role the most times.
ax = net.line_geodata.plot(column=net.line.cause_number, legend=True)
ax.set_title("Occurrence of causing the highest loading of other lines")
ax.set_axis_off();
The lines that cause the highest congestion in case of their failure can be represented graphically as demonstrated below. Also in this case the line 5 is the most critical.
ax = net.line_geodata.plot(column=net.res_line.congestion_caused_mva, legend=True)
ax.set_title("Congestion (MVA)")
ax.set_axis_off();