Tracing Infeasibilities¶

This tutorial demonstrates how to identify and trace infeasibilities in PyPSA optimization models using built-in functionality. When an optimization problem becomes infeasible, PyPSA provides tools to help you understand what constraints are conflicting and causing the infeasibility.

In this example, we'll deliberately create an infeasible network and then use PyPSA's built-in methods to diagnose the problem.

Getting ready¶

First, let's import PyPSA and load a example network from our example suite that we'll use for demonstration.

We'll use PyPSA's built-in AC-DC meshed example network, which includes:

• AC transmission lines connecting different regions
• A DC link providing additional transmission capacity
• Generators and loads across multiple buses

In [2]:

import pypsa

n = pypsa.examples.ac_dc_meshed()

n

import pypsa n = pypsa.examples.ac_dc_meshed() n

INFO:pypsa.network.io:Imported network 'AC-DC-Meshed' has buses, carriers, generators, global_constraints, lines, links, loads

Out[2]:

PyPSA Network 'AC-DC-Meshed'
----------------------------
Components:
 - Bus: 9
 - Carrier: 6
 - Generator: 6
 - GlobalConstraint: 1
 - Line: 7
 - Link: 4
 - Load: 6
Snapshots: 10

For this demonstration, we'll limit the analysis to just the first snapshot to keep the example simple and focused on the infeasibility tracing concept.

In [3]:

## Solve only the first period
n.snapshots = n.snapshots[:1]

## Solve only the first period n.snapshots = n.snapshots[:1]

This network is normally feasible, but let's modify it to create infeasibilities:

In [4]:

# Remove AC transmission lines connecting to London
n.remove("Line", "0")  # First AC line to London
n.remove("Line", "5")  # Second AC line to London

# Disable the DC link by setting capacity to zero
n.links.loc["DC link", "p_nom"] = 0.0
n.links.loc["DC link", "p_nom_extendable"] = False

# Remove AC transmission lines connecting to London n.remove("Line", "0") # First AC line to London n.remove("Line", "5") # Second AC line to London # Disable the DC link by setting capacity to zero n.links.loc["DC link", "p_nom"] = 0.0 n.links.loc["DC link", "p_nom_extendable"] = False

An alternative approach would be to remove the DC link entirely, which would also create an infeasible situation (London with load but no connections):

n.remove("Link", "DC link")

This would completely remove the DC link from the network. However, when n.optimize() is called, PyPSA would detect an "empty LHS and non-empty RHS" error during the constraint creation phase (specifically when building nodal balance constraints) and raise a ValueError before the optimization model is even passed to the solver.

Let's ensure that London has power load but no generators or connections:

In [5]:

print(f"London load: {n.loads_t.p_set.loc['2015-01-01', 'London']:.1f} MW")
print(
    f"London generation: {n.generators[n.generators.bus == 'London']['p_nom'].sum():.1f} MW"
)
print(
    f"London AC lines connected: {len(n.lines[(n.lines.bus0 == 'London') | (n.lines.bus1 == 'London')])}"
)
print(f"London DC link capacity: {n.links.loc['DC link', 'p_nom']:.1f} MW")
print("→ Infeasible: London has power load (RHS) but no supply options!")

print(f"London load: {n.loads_t.p_set.loc['2015-01-01', 'London']:.1f} MW") print( f"London generation: {n.generators[n.generators.bus == 'London']['p_nom'].sum():.1f} MW" ) print( f"London AC lines connected: {len(n.lines[(n.lines.bus0 == 'London') | (n.lines.bus1 == 'London')])}" ) print(f"London DC link capacity: {n.links.loc['DC link', 'p_nom']:.1f} MW") print("→ Infeasible: London has power load (RHS) but no supply options!")

London load: 35.8 MW
London generation: 0.0 MW
London AC lines connected: 0
London DC link capacity: 0.0 MW
→ Infeasible: London has power load (RHS) but no supply options!

Attempting Optimization¶

Now let's try to optimize this infeasible network. The optimization will fail because London has load but no way to meet it (no local generation and no transmission connections)

In [6]:

n.optimize(solver_name="gurobi")

n.optimize(solver_name="gurobi")

/tmp/ipykernel_8389/1458849580.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize(solver_name="gurobi")
WARNING:pypsa.consistency:The following lines have zero x, which could break the linear load flow:
Index(['2', '3', '4'], dtype='object', name='name')

WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['1', '6'], dtype='object', name='name')

INFO:linopy.model: Solve problem using Gurobi solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.07s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-hz8qvltj.lp

Reading time = 0.00 seconds

obj: 55 rows, 30 columns, 96 nonzeros

Set parameter LogToConsole to value 0

Warning: environment still referenced so free is deferred (Continue to use WLS)

WARNING:linopy.constants:Optimization potentially failed: 
Status: warning
Termination condition: infeasible_or_unbounded
Solution: 0 primals, 0 duals
Objective: nan
Solver model: available
Solver message: 4

Out[6]:

('warning', 'infeasible_or_unbounded')

Tracing the Infeasibility¶

When the optimization fails due to infeasibility, PyPSA provides a convenient method to diagnose the problem. The print_infeasibilities() method on the optimization model will show us exactly which constraints are causing the infeasibility.

In [7]:

n.model.print_infeasibilities()

n.model.print_infeasibilities()

Link-fix-p-lower[2015-01-01 00:00:00, DC link]: +1 Link-p[2015-01-01 00:00:00, DC link] ≥ -0
Bus-nodal_balance[London, 2015-01-01 00:00:00]: -1 Link-p[2015-01-01 00:00:00, DC link] = 35.7962441027

Interpreting the Infeasibility Output¶

The infeasibility trace above shows us exactly what's wrong:

1. Link-fix-p-lower[2015-01-01 00:00:00, DC link]: +1 Link-p[2015-01-01 00:00:00, DC link] ≥ -0

   • This constraint says the DC link power flow must be ≥ 0 (can't flow negative)
   • But we set the DC link capacity to 0, so Link-p value is bounded to 0

2. Bus-nodal_balance[London, 2015-01-01 00:00:00]: -1 Link-p[2015-01-01 00:00:00, DC link] = 35.7962441027

   • This is London's nodal balance equation: generation ± trade ± storage = load
   • London has no generation, 35.8 MW load, and no AC interconnectors
   • The DC link (the only remaining connection) can't provide any power because its capacity is 0

Under the hood, print_infeasibilities() uses Gurobi's infeasibility analysis to identify conflicting constraints. Irreducible Inconsistent Set (IIS) identified by Gurobi shows these two constraints are fundamentally incompatible:

• London needs power inflow via the DC link to meet its load
• But the DC link is constrained to zero power
• This creates an impossible situation that no solution can satisfy