Small correction in heat consumer for mode MF_TR

[dlohmeier]

the return temperature of the heat consumer is not fixed for the solver anymore, but calculated based on an exact estimate of qext from inlet and controlled return temperature

[JonasPfeiffer123]

Dear @dlohmeier,

I'm currently exploring the new return temperature control feature for heat_consumer, introduced in v0.11. Previously, I utilized a custom controller to adjust the controlled_mass_flow in heat_consumer based on a given output temperature. Now, I’m transitioning to the new approach, but my calculations fail under certain configurations.

I examined the provided functionality and corresponding tests in src/pandapipes/test/api/test_components/test_heat_consumer.py. While experimenting, I observed the following:

1. The original example works as expected, where one heat consumer is assigned treturn_k and the other controlled_mdot_kg_per_s.
2. Both heat consumers specifying controlled_mdot_kg_per_s (as in previous versions) also works fine.
3. However, when both heat consumers are assigned only treturn_k, the bidirectional calculation fails. Using the sequential mode also does not work. For clarity, qext_w is defined for both heat consumers.

Below is a minimal example demonstrating the issue:

def test_heat_consumer_result_extraction
    net = pp.create_empty_network("net", add_stdtypes=False, fluid="water")

    juncs = pp.create_junctions(
        net,
        nr_junctions=6,
        pn_bar=5,
        tfluid_k=286,
        system=["flow"] * 3 + ["return"] * 3,
        geodata=[
            (0, 0),    # Junction 0 (Startpunkt)
            (10, 0),   # Junction 1 (Vorlauf)
            (20, 0),   # Junction 2 (Vorlauf)
            (20, -10), # Junction 3 (Rücklauf)
            (10, -10), # Junction 4 (Rücklauf)
            (0, -10),  # Junction 5 (Endpunkt Rücklauf)
        ],
        name=[
            "Junction 0", # Junction 0
            "Junction 1", # Junction 1
            "Junction 2", # Junction 2
            "Junction 3", # Junction 3
            "Junction 4", # Junction 4
            "Junction 5" # Junction 5
        ]
    )
    pp.create_pipes_from_parameters(net, juncs[[0, 1, 3, 4]], juncs[[1, 2, 4, 5]], k_mm=0.1, length_km=1,
                                            diameter_m=0.1022, system=["flow"] * 2 + ["return"] * 2, alpha_w_per_m2k=0.4,
                                            text_k=273.15, name=[
            "Flow Pipe 1",  # Pipe von Source zu Flow Node 1
            "Flow Pipe 2",  # Pipe von Flow Node 1 zu Flow Node 2
            "Return Pipe 1", # Pipe von Return Node 1 zu Return Node 2
            "Return Pipe 2"  # Pipe von Return Node 2 zu Sink
        ])
    pp.create_circ_pump_const_pressure(net, juncs[-1], juncs[0], 5, 2, 85+273.15, type='auto', name="Pump 1")
    pp.create_heat_consumer(net, juncs[1], juncs[4], treturn_k=60+273.15, qext_w=7500, name="Heat Consumer 1")
    pp.create_heat_consumer(net, juncs[2], juncs[3], treturn_k=75+273.15, qext_w=7500, name="Heat Consumer 2")

    #hydraulics only to check for lookup heat transfer error
    pp.pipeflow(net, iter=3)

    pp.pipeflow(net, mode="bidirectional", iter=50) # played around with iter

    return net

Observations:

When both HeatConsumers are set to treturn_k with a defined qext_w, the bidirectional calculation fails. From my understanding of the new implementation, it should be possible to define output temperatures for all HeatConsumers without needing to set controlled_mdot_kg_per_s. Is my assumption correct? Or is at least one heat consumer with a set mass flow (controlled_mdot_kg_per_s) required?

Thank you in advance for your clarification and any advice you can provide!

[SimonRubenDrauz]

Hi Jonas,

I looked at your example and got it working. Two things I encountered here:
1.) You need to be careful with the initialisation of the junction temperature. In the first iteration this is the starting point. If your initial flow junction temperature is below your controlled temperature, you will face a problem. It is best to initialise your junction temperature with the flow temperature at the circ pump.
2.) The u-value of your pipes is quite high. The difference between the flow and return temperature in the second heat consumer is less than 3K. Due to the effect of a Newton raphson to oscillate, this could lead to the same effect as in 1 and cause your simulation to fail. The alpha value I used is a damping factor helping to reduce the oscillating effect of the NR.

Here is the code I used to make it work:

def test_heat_consumer_result_extraction():
net = pp.create_empty_network("net", add_stdtypes=False, fluid="water")

juncs = pp.create_junctions(
    net,
    nr_junctions=6,
    pn_bar=5,
    tfluid_k=85+273.15,
    system=["flow"] * 3 + ["return"] * 3,
    geodata=[
        (0, 0),    # Junction 0 (Startpunkt)
        (10, 0),   # Junction 1 (Vorlauf)
        (20, 0),   # Junction 2 (Vorlauf)
        (20, -10), # Junction 3 (Rücklauf)
        (10, -10), # Junction 4 (Rücklauf)
        (0, -10),  # Junction 5 (Endpunkt Rücklauf)
    ],
    name=[
        "Junction 0", # Junction 0
        "Junction 1", # Junction 1
        "Junction 2", # Junction 2
        "Junction 3", # Junction 3
        "Junction 4", # Junction 4
        "Junction 5" # Junction 5
    ]
)
pp.create_pipes_from_parameters(net, juncs[[0, 1, 3, 4]], juncs[[1, 2, 4, 5]], k_mm=0.1, length_km=1,
                                        diameter_m=0.1022, system=["flow"] * 2 + ["return"] * 2, u_w_per_m2k=0.4,
                                        text_k=273.15, name=[
        "Flow Pipe 1",  # Pipe von Source zu Flow Node 1
        "Flow Pipe 2",  # Pipe von Flow Node 1 zu Flow Node 2
        "Return Pipe 1", # Pipe von Return Node 1 zu Return Node 2
        "Return Pipe 2"  # Pipe von Return Node 2 zu Sink
    ])
pp.create_circ_pump_const_pressure(net, juncs[-1], juncs[0], 5, 2, 85+273.15, type='auto', name="Pump 1")
pp.create_heat_consumer(net, juncs[1], juncs[4], treturn_k=60+273.15, qext_w=7500, name="Heat Consumer 1")
pp.create_heat_consumer(net, juncs[2], juncs[3], treturn_k=75+273.15, qext_w=7500, name="Heat Consumer 2")

#hydraulics only to check for lookup heat transfer error
pp.pipeflow(net, iter=3)

pp.pipeflow(net, mode="bidirectional", iter=100, alpha=0.2) # played around with iter

return net

[JonasPfeiffer123]

Hi Simon,

Thank you very much for the quick help! The solution works perfectly for the simple test case.

Building on this, I encountered a similar issue during time-series calculations, which might stem from the same underlying problem, but I’m not entirely sure how to handle it in this context.

After verifying the simple network works, I returned to my more advanced network functions. The initialization and pipeflow calculation work flawlessly. However, during time-series testing—where I dynamically adjust qext_W at every time step—I ran into a problem. For most of the time steps, everything works fine, but in my dataset (generated using standard load profiles), there are some time steps where all heat_consumers have qext_W=0.

Previously, I handled this by setting a small minimum mass flow through the heat_consumer. While this workaround prevented the calculation from failing, it effectively ignored the desired output temperature and resulted in an artificially increased return temperature.

To illustrate the issue, here’s a chart showing this behavior (the chart is in German, but I think it’s clear enough). Heat demand and mass flow drop to zero, and the output temperature falls to 20 °C, which seems to be a default value when the heat_consumer is inactive (that 20 °C is the return joint to the circ_pump_const_pressure, which indicates that the whole net is 20 °C on the return line). However, in subsequent time steps, when qext_W rises again, the calculation fails to converge. This seems to be caused by the lower temperatures at the junctions, similar to the issue I experienced earlier.

I suspect that this default behavior (20 °C as a fallback temperature) may be the root cause here, but I couldn’t find the relevant part in the source code to confirm or modify this behavior.

Do you have any suggestions for handling this corner case? For example:

1. Is there a better way to initialize junction temperatures during these zero-demand time steps?
2. Could I enforce a minimum temperature at the junctions or some other constraint to stabilize the solver when all qext_W values are zero?
3. Are there alternative approaches for zero-demand time steps?

Thank you again for your support!

[SimonRubenDrauz]

Hi Jonas,

seems for me like a bug. However, without a code example, it is difficult to reproduce. Could you please provide an example.

[JonasPfeiffer123]

Hi Simon,

I'd like to, but currently the net where the problem occurs is way to data dependent (based on geojson input and custom net creation functions). I'm trying to get an easier example (not) working, where the same issue occurs. While doing that, I can show an example for net data.

My timeseries call looks like this:

def thermohydraulic_time_series_net(net, yearly_time_steps, qext_w_profiles, start, end, supply_temperature=85, supply_temperature_heat_consumer=75, return_temperature_heat_consumer=60):
    """Run a thermohydraulic time series simulation for the network.

    Args:
        net (pandapipesNet): The pandapipes network.
        yearly_time_steps (array): Array of yearly time steps.
        qext_w_profiles (list of arrays): List of external heat profiles.
        start (int): Start index for the simulation.
        end (int): End index for the simulation.
        supply_temperature (float, optional): Supply temperature. Defaults to 85.
        supply_temperature_heat_consumer (float, optional): Minimum supply temperature for heat consumers. Defaults to 75.
        return_temperature_heat_consumer (float, optional): Return temperature for heat consumers. Defaults to 60.

    Returns:
        tuple: Updated yearly time steps, network, and results.
    """
    # Prepare time series calculation
    yearly_time_steps = yearly_time_steps[start:end]

    # Update the ConstControl
    print(f"qext_w_profiles: {qext_w_profiles}") # Structure is two-dimensional array with shape (n_profiles, n_time_steps)
    # print minimum and maximum values of qext_w_profiles for each profile
    for i, qext_w_profile in enumerate(qext_w_profiles):
        print(f"qext_w_profile {i}: min={np.min(qext_w_profile)}, max={np.max(qext_w_profile)}")
    time_steps = range(0, len(qext_w_profiles[0][start:end]))
    update_const_controls(net, qext_w_profiles, time_steps, start, end)

    # If supply_temperature data exists, update corresponding Controller
    if supply_temperature is not None and isinstance(supply_temperature, np.ndarray):
        update_supply_temperature_controls(net, supply_temperature, time_steps, start, end)

    # Log variables and run time series calculation
    log_variables = create_log_variables(net)
    ow = OutputWriter(net, time_steps, output_path=None, log_variables=log_variables)

    run_time_series.run_timeseries(net, time_steps, mode="bidirectional", iter=100)

    return yearly_time_steps, net, ow.np_results

I do not use alpha=0.2 here, because when doing that I get the pandapipes.pf.pipeflow_setup.PipeflowNotConverged: The number of infeeding nodes and slacks do not match error in every time step.

In the first txt-file the results after a timeseries calculation are shown where all qext_w are >0 and no issue occurs.
net_results_working_time_step.txt

In the second file the results are shown for one time step later where some of the heat consumers have qext_w=0 which results in the calculation of NaN-Values in all connected junctions, but the calculation finishes.
net_results_not_working_time_step.txt

The problem starts to come up on the following time step when qext_w values would be higher again. I get the following output

dgstrf info 4
INFO:pandapipes.pf.pipeflow_setup:Setting the following nodes out of service for heat_transfer calculation in connectivity check:
In table junction: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87]
In table pipe_nodes: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343]

Which results in pandapipes.pf.pipeflow_setup.PipeflowNotConverged: The bidrectional calculation did not converge to a solution.

It's not really surprising that it fails, as half of the values are NaN. In the other time step I described in the previous comment where all qext_w are zero, it looks basically the same, only that all values are NaN.

Interestingly, when I test a small pipeflow calculation with one qext_w=0 value, which would look like the following, the values are also NaN for the heat consumer, but not NaN for the connected junctions. Also the Warning WARNING:pandapipes.component_models.heat_consumer_component:qext_w is equals to zero for heat consumers with index Index([1], dtype='int64'). Therefore, the defined temperature control cannot be maintained. occurs, which doesn't occur during the timeseries calculation:

      p_bar         t_k
0  2.500000  326.011735
1  4.000000  358.150000
2  3.985021  358.032986
3  3.910124  357.450604
4  3.872672  357.161085
5  2.973100  328.149999
6  2.591666  326.417565
7  2.515281  326.079110
   v_mean_m_per_s  p_from_bar  p_to_bar    t_from_k      t_to_k  t_outlet_k  \
0        0.459323    4.000000  3.985021  358.150000  358.032986  358.032986
1        0.459216    3.985021  3.910124  358.032986  357.450604  357.450604
2        0.459085    3.910124  3.872672  357.450604  357.161085  357.161085
3        0.451181    2.973100  2.591666  328.149999  326.417565  326.417565
4        0.450959    2.591666  2.515281  326.417565  326.079110  326.079110
5        0.450916    2.515281  2.500000  326.079110  326.011735  326.011735

   mdot_from_kg_per_s  mdot_to_kg_per_s  vdot_m3_per_s      reynolds    lambda
0            0.164532         -0.164532       0.000170  29019.235747  0.031813
1            0.164532         -0.164532       0.000170  28896.422165  0.031822
2            0.164532         -0.164532       0.000170  28738.919238  0.031834
3            0.164532         -0.164532       0.000167  18922.774136  0.032990
4            0.164532         -0.164532       0.000167  18598.550185  0.033048
5            0.164532         -0.164532       0.000167  18536.249740  0.033060
   p_from_bar  p_to_bar    t_from_k      t_to_k  t_outlet_k  \
0    3.872672    2.9731  357.161085  328.149999      328.15
1         NaN       NaN         NaN         NaN         NaN

   mdot_from_kg_per_s  mdot_to_kg_per_s  vdot_m3_per_s  reynolds  lambda  \
0            0.164532         -0.164532       0.000168       NaN     NaN
1                 NaN               NaN            NaN       NaN     NaN

    deltat_k   qext_w
0  29.011085  20000.0
1  64.300604      0.0
   v_mean_m_per_s  p_from_bar  p_to_bar    t_from_k  t_to_k  t_outlet_k  \
0             NaN         2.5       4.0  326.011735  358.15  326.011735

   mdot_from_kg_per_s  mdot_to_kg_per_s  vdot_m3_per_s  deltat_k  qext_w
0            0.164532         -0.164532       0.000167       NaN     NaN

So in conclusion the question would be why there are NaN values in the junction results in the timeseries calculation, when there are heat consumers with qext_w=0, which should result in an inactive heat consumer in that time step.

Maybe this description might give you an idea at what I'm looking at, I'll try to get an easy to test example running.

[JonasPfeiffer123]

Got a "minimal" example running where the same issue occurs. The following code includes all helper functions and classes I'm currently using. In the test function 20 time steps are calculated. When setting one qext_value=0 in the first heat consumer (here in time step 8) The calculation fails in time step 9 with

pandapower\timeseries\run_time_series.py", line 82, in pf_not_converged
    raise ts_variables['errors'][0] pandapipes.pf.pipeflow_setup.PipeflowNotConverged

The more important Information is coming from the output after time step 8:

"WARNING:pandapipes.component_models.heat_consumer_component:qext_w is equals to zero for heat consumers with index Index([0], dtype='int64'). Therefore, the defined temperature control cannot be maintained.
INFO:pandapipes.pf.pipeflow_setup:Setting the following nodes out of service for heat_transfer calculation in connectivity check:
In table junction: [1, 2, 3, 4, 5, 6, 7]
dgstrf info 8
C:\Users\jonas\AppData\Local\Programs\Python\Python311\Lib\site-packages\pandapipes\pipeflow.py:255: MatrixRankWarning: Matrix is exactly singular
  x = spsolve(jacobian, epsilon)
INFO:pandapipes.pf.pipeflow_setup:Setting the following nodes out of service for heat_transfer calculation in connectivity check:
In table junction: [1, 2, 3, 4, 5, 6, 7]
dgstrf info 8"

In my understanding the nodes connected to the first heat_consumer are set out of service as the heat_consumer results in NaN-Values and as this are connection nodes to all other nodes this sets all other nodes out of service?

What am I doing wrong?

"""
Filename: 07_example_timeseries_pandapipes.py
Author: Dipl.-Ing. (FH) Jonas Pfeiffer
Date: 2024-12-06
Description: Script for testing the pandapipes net simulation functions.
Usage: Run the script to generate a simple pandapipes network.
Functions:
    initialize_test_net(qext_w=np.array([100000, 100000, 100000, 100000]),
                        return_temperature=np.array([55, 60, 50, 60]),
                        supply_temperature=85,
                        flow_pressure_pump=4,
                        lift_pressure_pump=1.5,
                        pipetype="110/202 PLUS")
    timeseries_test(net)
    create_controllers(net, qext_w)
    update_const_controls(net, qext_w_profiles, time_steps, start, end)
    create_log_variables(net)
Classes:
    WorstPointPressureController(BasicCtrl)
Example:
    $ python 07_example_timeseries_pandapipes.py
"""

import traceback
import logging
# Initialize logging
logging.basicConfig(level=logging.INFO)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import pandapipes as pp
import pandapipes.plotting as pp_plot
from pandapipes.timeseries import run_time_series
from pandapower.timeseries import OutputWriter
from pandapower.timeseries import DFData
from pandapower.control.controller.const_control import ConstControl
from pandapower.control.basic_controller import BasicCtrl

# from districtheatingsim.net_simulation_pandapipes.pp_net_time_series_simulation import update_const_controls, create_log_variables
# districtheatingsim.net_simulation_pandapipes.utilities import create_controllers

class WorstPointPressureController(BasicCtrl):
    """
    A controller for maintaining the pressure difference at the worst point in the network.
    
    Args:
        net (pandapipesNet): The pandapipes network.
        circ_pump_pressure_idx (int, optional): Index of the circulation pump. Defaults to 0.
        target_dp_min_bar (float, optional): Target minimum pressure difference in bar. Defaults to 1.
        tolerance (float, optional): Tolerance for pressure difference. Defaults to 0.2.
        proportional_gain (float, optional): Proportional gain for the controller. Defaults to 0.2.
        **kwargs: Additional keyword arguments.
    """
    def __init__(self, net, circ_pump_pressure_idx=0, target_dp_min_bar=1, tolerance=0.2, proportional_gain=0.2, **kwargs):
        super(WorstPointPressureController, self).__init__(net, **kwargs)
        self.circ_pump_pressure_idx = circ_pump_pressure_idx
        self.target_dp_min_bar = target_dp_min_bar
        self.tolerance = tolerance
        self.proportional_gain = proportional_gain

        self.iteration = 0  # Add iteration counter

    def calculate_worst_point(self, net):
        """Calculate the worst point in the heating network, defined as the heat exchanger with the lowest pressure difference.

        Args:
            net (pandapipesNet): The pandapipes network.

        Returns:
            tuple: The minimum pressure difference and the index of the worst point.
        """
        
        dp = []

        for idx, qext, p_from, p_to in zip(net.heat_consumer.index, net.heat_consumer["qext_w"], net.res_heat_consumer["p_from_bar"], net.res_heat_consumer["p_to_bar"]):
            if qext != 0:
                dp_diff = p_from - p_to
                dp.append((dp_diff, idx))

        # Find the minimum delta p where the heat flow is not zero
        dp_min, idx_min = min(dp, key=lambda x: x[0])

        return dp_min, idx_min

    def time_step(self, net, time_step):
        """Reset the iteration counter at the start of each time step.

        Args:
            net (pandapipesNet): The pandapipes network.
            time_step (int): The current time step.

        Returns:
            int: The current time step.
        """
        self.iteration = 0  # reset iteration counter
        self.dp_min, self.heat_consumer_idx = self.calculate_worst_point(net)

        return time_step

    def is_converged(self, net):
        """Check if the controller has converged.

        Args:
            net (pandapipesNet): The pandapipes network.

        Returns:
            bool: True if converged, False otherwise.
        """
        
        current_dp_bar = net.res_heat_consumer["p_from_bar"].at[self.heat_consumer_idx] - net.res_heat_consumer["p_to_bar"].at[self.heat_consumer_idx]

        # Check if the pressure difference is within tolerance
        dp_within_tolerance = abs(current_dp_bar - self.target_dp_min_bar) < self.tolerance

        if dp_within_tolerance == True:
            return dp_within_tolerance

    def control_step(self, net):
        """Adjust the pump pressure to maintain the target pressure difference.

        Args:
            net (pandapipesNet): The pandapipes network.
        """
        # Increment iteration counter
        self.iteration += 1

        # Check whether the heat flow in the heat exchanger is zero
        current_dp_bar = net.res_heat_consumer["p_from_bar"].at[self.heat_consumer_idx] - net.res_heat_consumer["p_to_bar"].at[self.heat_consumer_idx]
        current_plift_bar = net.circ_pump_pressure["plift_bar"].at[self.circ_pump_pressure_idx]
        current_pflow_bar = net.circ_pump_pressure["p_flow_bar"].at[self.circ_pump_pressure_idx]

        dp_error = self.target_dp_min_bar - current_dp_bar
        
        plift_adjustment = dp_error * self.proportional_gain
        pflow_adjustment = dp_error * self.proportional_gain        

        new_plift = current_plift_bar + plift_adjustment
        new_pflow = current_pflow_bar + pflow_adjustment
        
        net.circ_pump_pressure["plift_bar"].at[self.circ_pump_pressure_idx] = new_plift
        net.circ_pump_pressure["p_flow_bar"].at[self.circ_pump_pressure_idx] = new_pflow

        return super(WorstPointPressureController, self).control_step(net)

def create_controllers(net, qext_w):
    """Create controllers for the network to manage heat consumers.

    Args:
        net (pandapipesNet): The pandapipes network.
        qext_w (array-like): External heat values for heat consumers.

    Returns:
        pandapipesNet: The pandapipes network with controllers added.
    """

    # Creates controllers for the network
    for i in range(len(net.heat_consumer)):
        # Create a simple DFData object for qext_w with the specific value for this pass
        placeholder_df = pd.DataFrame({f'qext_w_{i}': [qext_w[i]]})
        placeholder_data_source = DFData(placeholder_df)

        ConstControl(net, element='heat_consumer', variable='qext_w', element_index=i, data_source=placeholder_data_source, profile_name=f'qext_w_{i}')

    #dp_min, idx_dp_min = calculate_worst_point(net)  # This function must be defined
    dp_controller = WorstPointPressureController(net)#, idx_dp_min)  # This class must be defined
    net.controller.loc[len(net.controller)] = [dp_controller, True, -1, -1, False, False]

    return net

def update_const_controls(net, qext_w_profiles, time_steps, start, end):
    """Update constant controls with new data sources for time series simulation.

    Args:
        net (pandapipesNet): The pandapipes network.
        qext_w_profiles (list of arrays): List of external heat profiles.
        time_steps (range): Range of time steps for the simulation.
        start (int): Start index for slicing the profiles.
        end (int): End index for slicing the profiles.
    """
    for i, qext_w_profile in enumerate(qext_w_profiles):
        df = pd.DataFrame(index=time_steps, data={f'qext_w_{i}': qext_w_profile[start:end]})
        data_source = DFData(df)
        for ctrl in net.controller.object.values:
            if isinstance(ctrl, ConstControl) and ctrl.element_index == i and ctrl.variable == 'qext_w':
                ctrl.data_source = data_source

def create_log_variables(net):
    """Create a list of variables to log during the time series simulation.

    Args:
        net (pandapipesNet): The pandapipes network.

    Returns:
        list: List of tuples representing the variables to log.
    """
    log_variables = [
        ('res_junction', 'p_bar'), 
        ('res_junction', 't_k'),
        ('heat_consumer', 'qext_w'),
        ('res_heat_consumer', 'vdot_m3_per_s'),
        ('res_heat_consumer', 't_from_k'),
        ('res_heat_consumer', 't_to_k'),
        ('res_heat_consumer', 'mdot_from_kg_per_s'),
        ('res_circ_pump_pressure', 'mdot_from_kg_per_s'),
        ('res_circ_pump_pressure', 'p_to_bar'),
        ('res_circ_pump_pressure', 'p_from_bar')
    ]

    if 'circ_pump_mass' in net:
        log_variables.append(('res_circ_pump_mass', 'mdot_from_kg_per_s'))
        log_variables.append(('res_circ_pump_mass', 'p_to_bar'))
        log_variables.append(('res_circ_pump_mass', 'p_from_bar'))

    return log_variables

def initialize_test_net(qext_w=np.array([100000, 100000, 100000]),
                        return_temperature=np.array([55, 60, 50]),
                        supply_temperature=85,
                        flow_pressure_pump=4, 
                        lift_pressure_pump=1.5,
                        pipetype="110/202 PLUS"):
    net = pp.create_empty_network(fluid="water")
    
    k = 0.1  # roughness
    supply_temperature_k = supply_temperature + 273.15  # convert to Kelvin
    return_temperature_k = return_temperature + 273.15  # convert to Kelvin

    # Define junctions
    j1 = pp.create_junction(net, pn_bar=1.05, tfluid_k=supply_temperature_k, name="Pump Supply", geodata=(0, 0))
    j2 = pp.create_junction(net, pn_bar=1.05, tfluid_k=supply_temperature_k, name="Main Split Supply", geodata=(10, 0))
    j12 = pp.create_junction(net, pn_bar=1.05, tfluid_k=supply_temperature_k, name="Main Split Return", geodata=(10, 10))
    j13 = pp.create_junction(net, pn_bar=1.05, tfluid_k=supply_temperature_k, name="Pump Return", geodata=(0, 10))

    # Additional junctions for new branches
    j3 = pp.create_junction(net, pn_bar=1.05, tfluid_k=supply_temperature_k, name="Consumer B Supply", geodata=(20, 0))
    j4 = pp.create_junction(net, pn_bar=1.05, tfluid_k=supply_temperature_k, name="Consumer B Return", geodata=(20, 10))
    j5 = pp.create_junction(net, pn_bar=1.05, tfluid_k=supply_temperature_k, name="Consumer C Supply", geodata=(30, 0))
    j6 = pp.create_junction(net, pn_bar=1.05, tfluid_k=supply_temperature_k, name="Consumer C Return", geodata=(30, 10))

    # Pump
    pp.create_circ_pump_const_pressure(net, j13, j1, p_flow_bar=flow_pressure_pump, plift_bar=lift_pressure_pump, 
                                       t_flow_k=supply_temperature_k, type="auto", name="Main Pump")

    # Pipes for supply line
    pp.create_pipe(net, j1, j2, std_type=pipetype, length_km=0.02, k_mm=k, name="Main Pipe Supply")
    pp.create_pipe(net, j2, j3, std_type=pipetype, length_km=0.03, k_mm=k, name="Branch B Pipe Supply")
    pp.create_pipe(net, j3, j5, std_type=pipetype, length_km=0.03, k_mm=k, name="Branch C Pipe Supply")

    # Pipes for return line
    pp.create_pipe(net, j12, j13, std_type=pipetype, length_km=0.02, k_mm=k, name="Main Pipe Return")
    pp.create_pipe(net, j4, j12, std_type=pipetype, length_km=0.03, k_mm=k, name="Branch B Pipe Return")
    pp.create_pipe(net, j6, j4, std_type=pipetype, length_km=0.03, k_mm=k, name="Branch C Pipe Return")

    # Heat consumers
    pp.create_heat_consumer(net, from_junction=j2, to_junction=j12, qext_w=qext_w[0], treturn_k=return_temperature_k[0], name="Consumer A")
    pp.create_heat_consumer(net, from_junction=j3, to_junction=j4, qext_w=qext_w[1], treturn_k=return_temperature_k[1], name="Consumer B")
    pp.create_heat_consumer(net, from_junction=j5, to_junction=j6, qext_w=qext_w[2], treturn_k=return_temperature_k[2], name="Consumer C")

    # Simulate pipe flow
    pp.pipeflow(net, mode="bidirectional", iter=100, alpha=0.2)

    # Placeholder functions for additional processing
    net = create_controllers(net, qext_w)

    return net


def timeseries_test(net):
    start = 0
    end = 20 # 8760 hours in a year

    # time steps with start and end
    time_steps = np.arange(start, end, 1)	# time steps in hours

    # yearly time steps with dates beginning at 01.01.2021 00:00:00
    yearly_time_steps = pd.date_range(start="2021-01-01 00:00:00", periods=end, freq="H")

    # np.random.seed() is used to make the random numbers predictable
    np.random.seed(0)
    # for every time step for every heat consumer qext_w needs to be defined and saved in a two-dimensional array, not zeros random numbers in range 0 to 100000
    qext_w_profiles = np.random.randint(0, 100000, size=(4, end)) # Structure is two-dimensional array with shape (n_profiles, n_time_steps)
    print(f"qext_w_profiles: {qext_w_profiles}") # Structure is two-dimensional array with shape (n_profiles, n_time_steps)

    # set some time steps to zero
    qext_w_profiles[0, 7] = 0
    #qext_w_profiles[1, 7] = 0
    #qext_w_profiles[2, 13] = 0

    #qext_w_profiles[2, 306] = 1000

    print(f"qext_w_profiles: {qext_w_profiles}") # Structure is two-dimensional array with shape (n_profiles, n_time_steps)

    print(f"qext_w_profiles[0][8]: {qext_w_profiles[0][8]}") # Structure is two-dimensional array with shape (n_profiles, n_time_steps)
    print(f"qext_w_profiles[1][8]: {qext_w_profiles[1][8]}") # Structure is two-dimensional array with shape (n_profiles, n_time_steps)
    print(f"qext_w_profiles[2][8]: {qext_w_profiles[2][8]}") # Structure is two-dimensional array with shape (n_profiles, n_time_steps)

    update_const_controls(net, qext_w_profiles, time_steps, start, end)

    # Log variables and run time series calculation
    log_variables = create_log_variables(net)
    ow = OutputWriter(net, time_steps, output_path=None, log_variables=log_variables)

    run_time_series.run_timeseries(net, time_steps, mode="bidirectional", iter=100, alpha=0.2)

    return yearly_time_steps, net, ow.np_results


if __name__ == "__main__":
    try:
        print("Running the heat consumer result extraction test.")
        net = initialize_test_net()

        print(net)
        print(net.junction)
        print(net.pipe)
        print(net.heat_consumer)
        print(net.circ_pump_pressure)

        print(net.res_junction)
        print(net.res_pipe)
        print(net.res_heat_consumer)
        print(net.res_circ_pump_pressure)

        # create ax for config_plot
        pp_plot.simple_plot(net, junction_size=0.01, heat_consumer_size=0.1, pump_size=0.1, 
                            pump_color='green', pipe_color='black', heat_consumer_color="blue", 
                            show_plot=False)

        yearly_time_steps, net, np_results = timeseries_test(net)

        #plt.plot()

    except Exception as e:
        print("An error occurred:")
        print(traceback.format_exc())
        raise e