Violation of energy conservation in case of dangling links

[ekatef]

Checklist

• I am using the current main branch or the latest release. Please indicate.
• I am running on an up-to-date pypsa-earth environment. Update via conda env update -f envs/environment.yaml.

Describe the Bug

Network topology issues can lead to the fact that the energy conservation is not being being maintained in the model. The bug was noticed in previously (e.g. #1501) and investigated in details in #1500, and a minimal reproducible example has been created by @hazemakhalek and shared in this comment. Copy-pasting it bellow for convenience with set_seed() added.

Short-term solution

Make sure that you don't have links in your model whose buses are not defined. The easiest way to check it is to load your network and make sure that there is no warnings like The following links have buses which are not defined.

Long-term solution

The issue is resolved in the current PyPSA version (v0.35.1) which throws an error in case there are dangling branches in the network.

Example

import pypsa
import random
import numpy as np

random.seed(123)
np.random.seed(123)

nt = pypsa.Network()
hours = 24  
nt.set_snapshots(range(24))

no_of_buses = 4  # Define the number of buses
link_p_min_pu = -1  # Minimum power flow as a fraction of p_nom

for i in range(no_of_buses):
    nt.add("Bus", f"Bus {i+1}")

# Connect buses with links 
for i in range(4):
    nt.add("Link",
           f"Link {i+1}",
           bus0=f"Bus {i+1}",
           bus1=f"Bus {(i+1)%5 + 1}",
           p_nom=100,
           p_min_pu=link_p_min_pu,)  # bidirectionality


# Add generators with random parameters to each bus
for i in range(4):
    nt.add("Generator",
           f"Gen {i+1}",
           bus=f"Bus {i+1}",
           p_nom=random.uniform(50, 150),  
           marginal_cost=random.uniform(20, 80), 
           efficiency=random.uniform(0.3, 0.6))

# Add different load profiles to each bus
for i in range(3):
    # Create a daily load profile with randomness
    base_load = random.uniform(30, 70)
    profile = base_load + 10 * np.sin(np.linspace(0, 2 * np.pi, hours)) + np.random.normal(0, 3, hours)
    nt.add("Load",
           f"Load {i+1}",
           bus=f"Bus {i+1}",
           p_set=profile)
    
nt.consistency_check()
nt.optimize()
#nt.lopf()