Error in .lp output leads to infeasibility

[koen-vg]

I have a model which is "infeasible" when the "lp" or "lp-polars" io APIs are used, but is feasible when using the "direct" API. I'm pretty sure that the "correct" answer is that the model is feasible. See this example:

>>> import linopy
>>> linopy.read_netcdf("model.nc")
Linopy LP model
===============

Variables:
----------
 * Generator-p (snapshot, Generator)
[...]

Constraints:
------------
 * Bus-meshed-nodal_balance (Bus-meshed, snapshot)
[...]

Status:
-------
warning
>>> m = linopy.read_netcdf("model.nc")
>>> m.solve()
Writing constraints.: 100%|███████████████████████████████████████| 46/46 [00:08<00:00,  5.42it/s]
Writing continuous variables.: 100%|██████████████████████████████| 15/15 [00:01<00:00, 12.48it/s]
Set parameter Username
Academic license - for non-commercial use only - expires 2024-12-04
Read LP format model from file /tmp/linopy-problem-n4ff15xx.lp
Reading time = 3.74 seconds
obj: 2362618 rows, 1076813 columns, 5593804 nonzeros
Gurobi Optimizer version 11.0.3 build v11.0.3rc0 (linux64 - "Fedora Linux 41 (Workstation Edition)")

CPU model: Intel(R) Core(TM) i9-9900K CPU @ 3.60GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 2362618 rows, 1076813 columns and 5593804 nonzeros
Model fingerprint: 0xb49ef40a
Coefficient statistics:
  Matrix range     [3e-03, 1e+03]
  Objective range  [2e-01, 2e+06]
  Bounds range     [9e-02, 4e+10]
  RHS range        [1e-01, 1e+09]
Warning: Model contains large rhs
Warning: Model contains large bounds
         Consider reformulating model or setting NumericFocus parameter
         to avoid numerical issues.
Presolve removed 91133 rows and 1344 columns
Presolve time: 0.39s

Solved in 0 iterations and 0.39 seconds (0.18 work units)
Infeasible or unbounded model
Optimization failed: 
Status: warning
Termination condition: infeasible_or_unbounded
Solution: 0 primals, 0 duals
Objective: nan
Solver model: available
Solver message: 4

('warning', 'infeasible_or_unbounded')



>>> m.solve(io_api="direct")
Set parameter Username
Academic license - for non-commercial use only - expires 2024-12-04
Gurobi Optimizer version 11.0.3 build v11.0.3rc0 (linux64 - "Fedora Linux 41 (Workstation Edition)")

CPU model: Intel(R) Core(TM) i9-9900K CPU @ 3.60GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 2309328 rows, 1076813 columns and 5593804 nonzeros
Model fingerprint: 0x6359e7ff
Coefficient statistics:
  Matrix range     [3e-03, 1e+03]
  Objective range  [2e-01, 2e+06]
  Bounds range     [9e-02, 4e+10]
  RHS range        [1e-01, 1e+09]
Warning: Model contains large rhs
Warning: Model contains large bounds
         Consider reformulating model or setting NumericFocus parameter
         to avoid numerical issues.
Presolve removed 1252290 rows and 131969 columns
Presolve time: 3.99s
Presolved: 1057038 rows, 944844 columns, 4121485 nonzeros

Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...

Elapsed ordering time = 10s
Ordering time: 21.60s

Barrier statistics:
 Dense cols : 1735
 Free vars  : 1871
 AA' NZ     : 7.230e+06
 Factor NZ  : 2.011e+08 (roughly 2.4 GB of memory)
 Factor Ops : 1.335e+12 (roughly 6 seconds per iteration)
 Threads    : 6

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   8.49113849e+15 -3.79704839e+15  1.27e+12 5.38e+04  6.03e+12    41s
   1   9.32023171e+15 -3.83222385e+15  1.14e+12 4.03e+04  4.84e+12    51s
   2   1.00621540e+16 -4.01490786e+15  1.02e+12 2.20e+04  3.50e+12    62s
   3   1.10278201e+16 -4.31533676e+15  7.28e+11 1.27e+04  2.35e+12    74s
[...] (solves to completion)

Here is a link to the model so you can replicate this: https://drive.google.com/file/d/1ih0Eh3KwJuDZko_zuADl88VXtzHyhF8Z/view?usp=sharing. Unfortunately it's on the large side (~630MB).

The model comes from pypsa-eur (master) with the following (very much default) configuration:

foresight: overnight

scenario:
  ll:
  - v1.5
  clusters:
  - 40
  sector_opts:
  - ""
  planning_horizons:
  - 2050

clustering:
  temporal:
    resolution_sector: "24h"

Unfortunately I cannot replicate the infeasibility on a smaller model (say, with only a subset of countries), it seems like the full spatial extent of pypsa-eur is needed. But the infeasibility is "robust" in pypsa-eur in that I really cannot get any serious instance of pypsa-eur to solve right now unless I set io_api="direct".

It took me a while to figure out that this might be a linopy problem, and while I now know it's linopy, I still don't have any idea how exactly this happens. Hopefully this can be figured out before too long so not too many pypsa-eur users have to struggle with inexplicable infeasibilities.

(I'm using linopy 4.0.0.)

[FabianHofmann]

mmh, not so good, I am wondering what is happening. In the meanwhile, I saw that the slice_size argument is not correctly passed down. could you do me a favor and check out #379 and set the slice_size in the solve function to None and try again?

[koen-vg]

Setting slice_size to None solves the problem! (Using #379 as instructed.) A shame I didn't get to test #365 since I was away :)

I wonder what kind of problem caused this since the regression apparently slipped through all the tests?

[FabianHofmann]

got it! the selection on the dimensions that are allowed to be sliced was not correctly processed. the slicing always happens on the largest dimension and in this case it mistakenly sliced the constraint on the term dimension which is the longest dimension.

[koen-vg]

Great! Glad this was solved in less time than it took me to triage the issue to the linopy io API :)

Might I suggest making a new point release with the fix rather immediately? As far as I can tell, any users of pypsa-eur that happen to install linopy=0.4.0 will run into problems.

[FabianHofmann]

closing this, v0.4.1 includes the bugfix :)