implement transformer thermal model

CHANGELOG.rst

@@ -10,6 +10,7 @@ Change Log
 - [CHANGED] ideal phase shifters only use ``tap_step_degree`` as angle per tap step. ``tap_step_percent`` must be NaN or 0
 - [REMOVED] ``tap_phase_shifter`` / ``tap2_phase_shifter``: phase shifters are only defined by ``tap_changer_type`` / ``tap2_changer_type``
 - [ADDED] added a sensitivity package, it provides LODF, PSDF and PTDF calculation based on the DC power flow
+- [ADDED] TDPF: top-oil and hot-spot temperatures of ``net.trafo`` and ``net.trafo3w``, based on IEC 60076-7 and the optional ``transformer-thermal-model`` package (``pip install pandapower[tdpf]``); the thermal model runs as post-processing on the solved loading, so ``r_sc`` is unchanged and voltages/flows do not change; new input columns ``cooling_type`` (ONAN/ONAF, default ONAN) and ``air_temperature_degree_celsius`` (default 35) and new result columns ``top_oil_temperature_degree_celsius``, ``hot_spot_temperature_degree_celsius`` (2W) / ``hot_spot_temperature_degree_celsius_hv|mv|lv`` (3W) and ``temperature_degree_celsius``; without the package a warning is logged and line TDPF is unaffected
 - [FIXED] ``res_dcline``: a terminal at an unsupplied bus (island without slack) now has ``NaN`` as voltage magnitude and angle, consistent with ``res_bus``, instead of ``0``
 - [FIXED] branch voltage angle difference constraints (MATPOWER ``ANGMIN``/``ANGMAX``) were lost in the pypower/matpower converters and ignored by the OPF: the new optional columns ``min_angle_degree`` and ``max_angle_degree`` in ``net.line``, ``net.trafo``, ``net.trafo3w`` and ``net.impedance`` are considered by ``runopp``/``rundcopp`` and by ``to_pm``, and ``from_ppc``/``to_ppc``/``to_mpc`` preserve them (#2385). For ``net.trafo3w`` the limit bounds each internal star branch (hv -> star and star -> mv/lv) separately, so the terminal difference can reach twice the configured value. OPF results of networks imported from MATPOWER/PYPOWER cases with angle limits change accordingly; pass ``OPF_IGNORE_ANG_LIM=True`` to ``runopp`` to ignore the limits.
 - [FIXED] matpower/pypower converter: transformers without a tap changer (tap ratio 1) were converted to impedance elements by a to_mpc/from_mpc (and to_ppc/from_ppc) round trip, since a tap ratio of 1 was not recognized as a transformer when importing. The converters now follow the matpower convention (a tap ratio of 0 denotes "no transformer", any other ratio denotes a transformer): ``to_mpc`` keeps the tap ratio of transformer branches, ``from_mpc`` no longer rewrites tap ratios of 0, ``to_ppc`` marks branches of impedance elements with a tap ratio of 0, and ``from_ppc`` only imports branches without a tap ratio and without a phase shift that connect different voltage levels as impedances. The tap ratio marking also stays correct when branches are filtered out of the internal ppc (out-of-service or disconnected branches). A transformer with a tap ratio of 1 between buses of the same voltage level still cannot be distinguished from a line in the ppc/matpower format and keeps being converted to a line on the way back, as before (#2620)

doc/about/license.rst

@@ -10,4 +10,3 @@ License
 pandapower is published under the following 3-clause BSD license:
 
     .. literalinclude:: ../../LICENSE
-

doc/elements/trafo.rst

@@ -35,6 +35,11 @@ Table Structure
 
 .. note:: The transformer loading constraint for the optimal power flow corresponds to the option trafo_loading="current":
 
+.. note:: The columns ``cooling_type`` and ``air_temperature_degree_celsius`` are only used by the
+    transformer thermal calculation of the temperature dependent power flow, see
+    :ref:`Transformer temperatures <tdpf_transformer_temperatures>`. They default to ``"ONAN"`` (Oil Natural, Air Natural) and
+    ``35`` °C, so they can be left empty.
+
 .. note:: vkr_percent can be calculated as follow :
 
     .. math::
@@ -404,6 +409,10 @@ Result Tables
     :header-rows: 1
     :delim: ,
 
+.. note:: With ``runpp(tdpf=True)`` the columns ``top_oil_temperature_degree_celsius``,
+    ``hot_spot_temperature_degree_celsius`` and ``temperature_degree_celsius`` are added, see
+    :ref:`Transformer temperatures <tdpf_transformer_temperatures>`.
+
 .. math::
    :nowrap:
 

doc/elements/trafo3w.rst

@@ -44,6 +44,11 @@ Table Structure
 
 .. note:: Three Winding Transformer loading can not yet be constrained with the optimal power flow.
 
+.. note:: The columns ``cooling_type`` and ``air_temperature_degree_celsius`` are only used by the
+    transformer thermal calculation of the temperature dependent power flow, see
+    :ref:`Transformer temperatures <tdpf_transformer_temperatures>`. They default to ``"ONAN"`` (Oil Natural, Air Natural) and
+    ``35`` °C, so they can be left empty.
+
 .. _trafo3w_model:
 
 Electric Model
@@ -173,6 +178,11 @@ Result Tables
     :header-rows: 1
     :delim: ,
 
+.. note:: With ``runpp(tdpf=True)`` the columns ``top_oil_temperature_degree_celsius``,
+    ``hot_spot_temperature_degree_celsius_hv``, ``hot_spot_temperature_degree_celsius_mv``,
+    ``hot_spot_temperature_degree_celsius_lv`` and ``temperature_degree_celsius`` are added, see
+    :ref:`Transformer temperatures <tdpf_transformer_temperatures>`.
+
 .. math::
     p\_hv\_mw &= Re(\underline{v}_{hv} \cdot \underline{i}_{hv}) \\    
     q\_hv\_mvar &= Im(\underline{v}_{hv} \cdot \underline{i}_{hv}) \\

doc/powerflow/ac.rst

@@ -127,6 +127,88 @@ More detailed information about TDPF can be found in the tutorials:
 
 - `TDPF introduction <https://github.com/e2nIEE/pandapower/blob/develop/tutorials/temperature_dependent_power_flow.ipynb>`__
 - `Sensitivity of TDPF parameters <https://github.com/e2nIEE/pandapower/blob/develop/tutorials/temperature_dependent_power_flow_parameters.ipynb>`__
+- `Transformer temperatures <https://github.com/e2nIEE/pandapower/blob/develop/tutorials/temperature_dependent_power_flow_trafo.ipynb>`__
+
+.. _tdpf_transformer_temperatures:
+
+Transformer temperatures
+````````````````````````
+
+In addition to the conductor temperature of the lines, TDPF reports the top-oil and hot-spot temperature
+of every in-service transformer (:code:`net.trafo` and :code:`net.trafo3w`). The calculation is based on
+the IEC 60076-7 loading guide and utilizes the python package
+`transformer-thermal-model <https://github.com/alliander-opensource/transformer-thermal-model>`__ developed by Alliander.
+
+This is an optional feature and requires an optional dependency. Install it with:
+
+.. code-block:: bash
+
+    pip install pandapower[tdpf]
+
+If the package is not installed, TDPF logs a warning and skips transformer temperature calculations.
+Line TDPF continues to work as before. Transformer temperature result columns are not created,
+so line-only TDPF users do not see columns filled with NaN values.
+
+Two additional input columns are read from :code:`net.trafo` and :code:`net.trafo3w`:
+
+.. tabularcolumns:: |p{0.30\linewidth}|p{0.15\linewidth}|p{0.45\linewidth}|
+.. csv-table::
+    :header: "Column", "Type", "Description"
+    :widths: 30, 15, 45
+
+    ``cooling_type``,string,"cooling type of the transformer, either ONAN or ONAF (default: ONAN)"
+    ``air_temperature_degree_celsius``,float,"ambient air temperature in °C (default: 35)"
+
+The temperature is *reported*, not fed back into the power flow. In contrast to the line TDPF, the
+transformer resistance :code:`r_sc` is never modified, because the thermal model runs as post-processing
+on the solved transformer loading. The results are therefore not iterated with the power flow and do not
+change voltages, flows or losses.
+
+The results are written to additional columns of :code:`net.res_trafo` and :code:`net.res_trafo3w`:
+
+.. tabularcolumns:: |p{0.35\linewidth}|p{0.30\linewidth}|p{0.35\linewidth}|
+.. csv-table::
+    :header: "Column", "Table", "Description"
+    :widths: 35, 30, 35
+
+    ``top_oil_temperature_degree_celsius``,"res_trafo, res_trafo3w","top-oil temperature in °C"
+    ``hot_spot_temperature_degree_celsius``,"res_trafo","hot-spot temperature in °C"
+    ``hot_spot_temperature_degree_celsius_hv``,"res_trafo3w","hot-spot temperature of the HV winding in °C"
+    ``hot_spot_temperature_degree_celsius_mv``,"res_trafo3w","hot-spot temperature of the MV winding in °C"
+    ``hot_spot_temperature_degree_celsius_lv``,"res_trafo3w","hot-spot temperature of the LV winding in °C"
+    ``temperature_degree_celsius``,"res_trafo, res_trafo3w","hottest hot-spot temperature in °C, used by the contingency limit checks"
+
+For a three winding transformer, :code:`temperature_degree_celsius` is the maximum of the three winding
+hot-spot temperatures, mirroring how :code:`loading_percent` aggregates the windings.
+
+Whether the model is solved for the steady state or for a transient depends on the :code:`tdpf_delay_s`
+power flow parameter, the same parameter that controls the thermal inertia of the lines:
+
+* :code:`tdpf_delay_s=None` (the default) gives the steady state temperature for the current loading.
+* A finite :code:`tdpf_delay_s` gives the temperature after this time step, starting from a cold
+  transformer. In a time series run (:code:`pandapower.timeseries.run_timeseries`), the top-oil temperature
+  of the previous time step is used as the initial condition, so the transformer warms up and cools down
+  continuously over the simulation.
+
+.. note::
+
+    The ONAN and ONAF cooling types share the same steady state temperature. They only differ in the
+    transient behaviour: an ONAF transformer has shorter oil and winding time constants and therefore
+    reaches its final temperature faster, not at a lower temperature.
+
+.. note::
+
+    Only the nameplate data that pandapower already stores is passed to the thermal model: the no-load
+    loss, the load loss, the rated apparent power and the rated voltages of each winding, which yield
+    the nominal currents. The optional specifications of the loading guide -- above all
+    :code:`winding_oil_gradient`, :code:`hot_spot_fac` and the winding and oil time constants -- are
+    left at the IEC 60076-7 defaults of the library (17 K, 1.3 and the values for the selected cooling
+    type).
+
+    The top-oil temperature depends only on the nameplate data and therefore matches a calculation with
+    modified optional specifications exactly. The hot-spot temperature does depend on them, so it will
+    differ from a calculation that uses a transformer-specific winding-oil gradient or hot-spot factor.
+    A wider set of input columns is not part of this implementation.
 
 .. autofunction:: pandapower.run.runpp
 

mypy.ini

@@ -64,5 +64,8 @@ ignore_missing_imports = true
 [mypy-pypower.*]
 ignore_missing_imports = true
 
+[mypy-transformer_thermal_model.*]
+ignore_missing_imports = true
+
 [mypy-xlsxwriter.*]
 ignore_missing_imports = true

pandapower/auxiliary.py

@@ -1500,6 +1500,10 @@
     if len(tdpf_lines) == 0:
         logger.info("TDPF: no relevant lines found")
 
+    # transformers are always included in the thermal calculation when tdpf=True (no opt-in flag)
+    trafo_idx = net.trafo.index[net.trafo.in_service.values] if len(net.trafo) else net.trafo.index[:0]
+    trafo3w_idx = net.trafo3w.index[net.trafo3w.in_service.values] if len(net.trafo3w) else net.trafo3w.index[:0]
+
     # required for simplified approach if r_theta_kelvin_per_mw is provided, can be filled with simplified values:
     default_values = {
         "temperature_degree_celsius": 20,  # starting temperature of lines
@@ -1533,20 +1537,22 @@
     # if np.any(np.setdiff1d(net.line.columns, required_columns)):
     #     raise UserWarning(f"TDPF: required columns missing in net.line: {required_columns}")
 
-    missing_columns = []
-    for col in required_columns:
-        if col not in net.line.columns or np.any(net.line.loc[tdpf_lines, col].isnull()):
-            missing_columns.append(col)
-
-    if len(missing_columns) > 0:
-        raise UserWarning(f"TDPF: required columns {missing_columns} are missing or have missing values")
-
-    # check if unsupported elements are included
-    if (
-        len(net.line.loc[net.line.index.isin(tdpf_lines) & (net.line.type == "cs")]) > 0
-        or "tdpf" in net.trafo
-        or "tdpf" in net.trafo3w
-    ):
+    # The line columns are only validated if the network contains lines at all. A network with only
+    # transformers (no lines) still gets transformer temperatures, so there is nothing to check here.
+    if len(net.line) > 0:
+        missing_columns = []
+        for col in required_columns:
+            # `tdpf` selects the lines and therefore has to exist for the branch conversion, but only
+            # the values of the selected lines are relevant. An ordinary line that did not opt in must
+            # not trigger a missing-value warning for a column it never needed.
+            if col not in net.line.columns or (len(tdpf_lines) > 0 and np.any(net.line.loc[tdpf_lines, col].isnull())):
+                missing_columns.append(col)
+
+        if len(missing_columns) > 0:
+            raise UserWarning(f"TDPF: required columns {missing_columns} are missing or have missing values")
+
+    # check if unsupported elements are included; the "type" column is optional in the line schema
+    if len(tdpf_lines) > 0 and "type" in net.line.columns and (net.line.loc[tdpf_lines, "type"].dropna() == "cs").any():
         logger.warning("TDPF: temperature dependent power flow is only implemented for overhead lines")
 
     # now fill in the default values
@@ -1557,9 +1563,46 @@
             logger.info(f"TDPF: filling nan values in {col} with a default assumption of {val}")
             net.line.loc[tdpf_lines, col] = net.line.loc[tdpf_lines, col].fillna(val)
 
-    if len(net.line.loc[net.line.index.isin(tdpf_lines) & (net.line.r_ohm_per_km == 0)]) > 0:
+    if len(net.line) > 0 and len(net.line.loc[net.line.index.isin(tdpf_lines) & (net.line.r_ohm_per_km == 0)]) > 0:
         raise UserWarning("TDPF: temperature dependent power flow cannot be applied to lines that have r_ohm_per_km=0")
 
+    _check_tdpf_trafo_parameters(net, trafo_idx, trafo3w_idx)
+
+
+def _check_tdpf_trafo_parameters(net: pandapowerNet, trafo_idx, trafo3w_idx) -> None:
+    """Default and validate the transformer input columns of the transformer thermal calculation.
+
+    Every in-service transformer takes part in the thermal calculation when ``tdpf=True``; there is no
+    per transformer opt-in flag. Only ``cooling_type`` and ``air_temperature_degree_celsius`` are
+    required inputs, everything else is derived from the nameplate data.
+    """
+    trafo_defaults = {"cooling_type": "ONAN", "air_temperature_degree_celsius": 35}
+    for element, indices in (("trafo", trafo_idx), ("trafo3w", trafo3w_idx)):
+        if len(indices) == 0:
+            continue
+        table = net[element]
+        if "cooling_type" not in table.columns:
+            table["cooling_type"] = pd.Series(pd.NA, index=table.index, dtype="string")
+        if "air_temperature_degree_celsius" not in table.columns:
+            table["air_temperature_degree_celsius"] = np.nan
+        for col, val in trafo_defaults.items():
+            if np.any(pd.isnull(table.loc[indices, col])):
+                logger.info(f"TDPF: filling nan values in {element}.{col} with a default assumption of {val}")
+        # normalise the cooling type, so that the library mapping is unambiguous
+        table["cooling_type"] = table["cooling_type"].astype("string")
+        table.loc[indices, "cooling_type"] = table.loc[indices, "cooling_type"].str.strip().str.upper().fillna("ONAN")
+        table.loc[indices, "air_temperature_degree_celsius"] = (
+            table.loc[indices, "air_temperature_degree_celsius"].astype(float).fillna(35.0)
+        )
+        invalid = ~table.loc[indices, "cooling_type"].isin(["ONAN", "ONAF"])
+        if invalid.any():
+            raise UserWarning(
+                f"TDPF: invalid cooling_type in {element} for indices {table.loc[indices].index[invalid].tolist()}, "
+                "must be 'ONAN' or 'ONAF'"
+            )
+        if (table.loc[indices, "air_temperature_degree_celsius"] <= -273).any():
+            raise UserWarning(f"TDPF: air_temperature_degree_celsius in {element} must be greater than -273 °C")
+
 
 # =============================================================================
 # Functions for 3 Phase Unbalanced Load Flow

pandapower/contingency/contingency.py

@@ -22,6 +22,7 @@
 
 from pandapower.create._utils import add_column_to_df
 from pandapower.network import pandapowerNet
+from pandapower.pf.transformer_thermal import TRANSFORMER_THERMAL_AVAILABLE
 from pandapower.run import runpp
 
 logger = logging.getLogger(__name__)
@@ -154,12 +155,23 @@ def run_contingency(
         "bus": ["vm_pu"],
         **{key: ["loading_percent"] for key in ("line", "trafo", "trafo3w") if len(net[key]) > 0},
     }
-    if len(net.line) > 0 and (
+    tdpf_requested = (
         net.get("_options", {}).get("tdpf", False)
         or pf_options.get("tdpf", False)
         or pf_options_nminus1.get("tdpf", False)
-    ):
-        result_variables["line"].append("temperature_degree_celsius")
+    )
+    if tdpf_requested:
+        if len(net.line) > 0:
+            result_variables["line"].append("temperature_degree_celsius")
+
+        # transformer temperatures are post-processing results of the transformer thermal model. They
+        # are produced by every solve once tdpf is on, no matter whether the input network already
+        # carries the result column -- a fresh network does not, so the column cannot be used to
+        # decide this before the first solve.
+        if TRANSFORMER_THERMAL_AVAILABLE:
+            for element in ("trafo", "trafo3w"):
+                if len(net[element]) > 0:
+                    result_variables[element].append("temperature_degree_celsius")
 
     # for n-1
     for element, val in nminus1_cases.items():
@@ -414,7 +426,11 @@ def _update_contingency_results(
 ):
     for element, vars in result_variables.items():
         for var in vars:
-            val = net[f"res_{element}"][var].values
+            res_table = net[f"res_{element}"]
+            if var not in res_table.columns:
+                # e.g. the transformer temperature when no transformer took part in this solve
+                continue
+            val = res_table[var].values
             if nminus1:
                 if var == "loading_percent":
                     s = (