Poor performance and NumericalError on gas pipeline network NLPs (PDE-discretized, banded structure)

[jkitchin]

Summary

Ripopt fails or runs very slowly on gas pipeline network NLP benchmark problems from Sakshi21299/gas_networks. These are Pyomo models built with the gas_net package on GasLib-11 and GasLib-40 test instances.

Two problem types arise from this package:

• Steady-state NLP: small, ~50–200 variables. Works fine.
• Dynamic NLP (24-hour horizon, finite-volume PDE discretization): n≈2500, m≈2500. Fails.

Benchmark Results (IPOPT vs RIPOPT)

Network    Phase    Solver   Status               Obj (kW·h)   Time (s)
------------------------------------------------------------------------
GasLib-11  steady   ipopt    optimal              0.0329       0.021
GasLib-11  steady   ripopt   optimal              0.0329       0.048   ✓

GasLib-11  dynamic  ipopt    optimal              0.7885       0.092
GasLib-11  dynamic  ripopt   internalSolverError  N/A          18.4    ✗

GasLib-40  steady   ipopt    optimal              1.2900       0.112
GasLib-40  steady   ripopt   internalSolverError  N/A          1246.6  ✗

Ripopt matches IPOPT exactly on the steady-state GasLib-11 problem. It fails on both the dynamic GasLib-11 and the steady-state GasLib-40.

Diagnostic Output (GasLib-11 dynamic, print_level=5)

ripopt: Starting main loop (n=2542, m=2492)
   0     1.7280000e0      8.64e2      0.00e0     -1.0    0.00e0    0.00e0    0
ripopt: Sparse condensed S has bandwidth 2541 for n=2542, switching to augmented KKT
   1     1.7281070e0      8.64e2      2.38e2     -2.0   1.81e-4   5.94e-4    0
ripopt: L-BFGS Hessian fallback did not improve (NumericalError)

--- ripopt diagnostics ---
status: NumericalError
iterations: 63
wall_time: 18.407s
final_mu: 1.01e-7
final_primal_inf: 2.21e-9
final_dual_inf: 2.96e-1
final_compl: 1.01e-7
---
ripopt 0.5.0: NumericalError after 63 iterations
Objective: 7.885488421941358e-1   ← matches IPOPT's optimal value exactly

Note: the final objective value matches IPOPT's optimal (0.7885), and final_primal_inf: 2.21e-9 is essentially zero. The solver found the right answer but dual convergence stalled.

Root Cause Hypothesis

The key line is:

Sparse condensed S has bandwidth 2541 for n=2542, switching to augmented KKT

The dynamic gas network NLP arises from a finite-volume discretization of the Euler equations on a pipe network. The Jacobian/Hessian has a block-banded structure (space × time grid), but ripopt's bandwidth heuristic is detecting bandwidth ≈ n, treating it as essentially dense and switching to the augmented KKT system.

Once in augmented KKT mode, the L-BFGS Hessian fallback is triggered and fails (NumericalError), causing dual convergence to stall even though the primal solution is correct.

The GasLib-40 steady-state case (n+m probably around 300–400) presumably crosses the sparse_threshold and hits the same bandwidth detection issue.

Problem Structure

The model is a Pyomo ConcreteModel built by gas_net/model_nlp.py. It uses:

• Finite-volume discretization of mass/momentum PDEs (backward time differencing)
• Epsilon-regularized flow reversal (|u| ≈ sqrt(u² + ε²))
• Scaled variables (pressure ÷ 1e5, power ÷ 1e5)
• Compressor performance maps as polynomial constraints

The .nl file produced by Pyomo likely has a variable ordering that makes the true banded structure non-obvious from the bandwidth of the condensed KKT matrix.

Options Tried

s.options['enable_lbfgs_hessian_fallback'] = 'no'   # prevents the fallback
s.options['warm_start_init_point'] = 'yes'
s.options['mu_init'] = 1e-6
s.options['bound_push'] = 1e-6

Disabling the L-BFGS fallback doesn't resolve the underlying dual convergence issue.

Suggested Directions

1. Variable reordering: The .nl file variable ordering may be scrambling the natural banded structure. A bandwidth-reducing reordering (AMD, RCM) applied before bandwidth detection might reveal the true structure.

2. Augmented KKT linear solver robustness: On the augmented system, the linear algebra appears to fail for this problem. Better inertia correction or a more robust factorization might help.

3. Dual convergence with near-zero primal infeasibility: When primal_inf < 1e-8 but dual_inf > 0.1, the solver is essentially at the primal solution but can't close the dual gap. An option to accept this as locallyOptimal (rather than NumericalError) with a warning would be useful — the current behavior discards a correct primal solution.

Reproduction

git clone https://github.com/Sakshi21299/gas_networks
pip install -e gas_networks/
pip install pyomo-ripopt

import json, pyomo.environ as pyo, pyomo_ripopt
from pathlib import Path
from gas_net.util.import_data import import_data_from_excel
from gas_net.model_nlp import buildNonLinearModel
from gas_net.modelling_library.fix_and_init_vars import init_network_default

DATA = Path("gas_networks/gas_net/data")
opt = json.load(open(DATA / "Options.json"))
opt["dynamic"] = True
opt["T"] = 24

nd, id_ = import_data_from_excel(
    str(DATA / "data_files/GasLib_11/networkData.xlsx"),
    str(DATA / "data_files/GasLib_11/inputData.xlsx"),
)
m = buildNonLinearModel(nd, id_, opt)
m = init_network_default(m, p_default=55e5)
m.wSource["entry01", :].fix()
m.wSource["entry02", :].fix()
m.pSource["entry03", :].fix()

s = pyo.SolverFactory("ripopt")
s.options["print_level"] = 5
result = s.solve(m, tee=True)

Environment

• ripopt 0.5.0 / 0.6.1
• Pyomo 6.x
• pyomo-ripopt 0.1.0
• macOS Darwin 25.3.0

[jkitchin]

Progress Update (2026-04-11)

Performance Fixes (gas problems now run)

Sparse GN restoration (src/restoration.rs): The Gauss-Newton restoration step was using dense Bunch-Kaufman factorization on J_a * J_a^T, which took ~6 seconds per step for m=2492. Replaced with sparse rmumps factorization for m>500, reducing to ~0.02s/step. This was the primary bottleneck preventing gaslib11_dynamic and gaslib40_steady from finishing.

KKT factorization is fast: rmumps factors the 5034x5034 augmented KKT in 0.009-0.21s. The bottleneck was never the KKT factorization itself.

Convergence Fix (Loqo mu oracle)

Root cause of the 164-iteration count (vs Ipopt's ~20): the mu oracle was disabled by default. The solver was using a simple fallback (avg_compl / kappa) that drops mu too aggressively from infeasible starting points, causing tiny fraction-to-boundary steps.

Ipopt keeps mu=0.1 for 10 iterations (its quality function oracle with golden-section finds sigma~1 when infeasible). Enabled the Loqo oracle (sigma = 0.1 * min(0.05*(1-xi)/xi, 2)^3) which is conservative when off-centered and aggressive when centered.

Results on gaslib11_steady:

Metric	Before	After	Ipopt
Iterations	164	134	20
NLP restorations	1	0	0
Mode switches	24	8	~3
Full steps (alpha=1)	rare	iters 7-9	iters 5-8

Remaining Gap (134 vs 20 iterations)

Based on Ipopt source analysis:

1. Quality function oracle: Ipopt does two extra KKT solves per iteration (affine + centering directions) and golden-section search over sigma. The Loqo oracle is a simpler approximation.

2. Linear solver quality: Ipopt/MUMPS uses MC64 scaling + pivtol=1e-6. rmumps uses Ruiz scaling + pivtol=0.01 (1e-6 regresses without MC64). Different search directions -> different convergence.

3. theta_min/theta_max: Ipopt sets once from initial theta; ripopt re-computes after each filter reset (workaround for slow convergence).

Status of gas problems

Problem	n+m	Before	After
gaslib11_steady	404	3.6s, 164 iters, Optimal	10.4s, 134 iters, Optimal
gaslib11_dynamic	5034	Never finished	Runs, IPM iters fast, NLP restoration slow
gaslib40_steady	3376	Never finished	Runs, IPM iters fast, NLP restoration slow

Next steps: full quality function oracle, rmumps MC64-equivalent scaling.