Extreme matrix range [6e-11, 2e+08] in multi-stage energy system MILP — model not solving within 24 hours

Hi there,

I am new to Gurobi and MILP and would be happy to get some input.

I am solving a multi-stage, multi-carrier energy system optimization model, solved with Gurobi 11.0. The model covers Switzerland's energy system across four stages (2020–2050) with ~1.9M rows, 750K columns, and 4.6M nonzeros.

Gurobi reports the following coefficient statistics:
Matrix range [6e-11, 2e+08]
Objective range [1e+00, 1e+00]
Bounds range [5e-03, 2e+08]
RHS range [7e-07, 1e+09]

I am aware this gives a matrix range ratio of ~3×10^18, far beyond Gurobi's recommended range. The model does not solve within 24 hours even with MIPGap=0.1 and only 110 TSAM-aggregated timesteps. Most of the time its stuck in the crossover phase.

I have few assumptions why matrix range is so big but I dont know how to fix it.

Solar PV capacity constraints include irradiance profiles in GW/m², where twilight TSAM clusters produce weighted coefficients as small as 9e-11. The system power unit is GW and area unit is m², making irradiance structurally small (~1e-7 GW/m² peak), but I do not know if changing units would help because other technologies like hydropower are just comparatively much more dominant in Switzerland than PV

I also have mixed mass/energy unit constraints (1e+06). For example, gasification technologies have both energy inputs (GWh) and mass inputs (tonnes of O2), producing input-ratio constraints with coefficients like -1,040,582 vs +25,641 in the same row.

Time series aggregation (TSAM) multiplies profile values by cluster weights (hours represented / total hours), further amplifying near-zero values.

Current solver settings*
NumericFocus = 2
ScaleFlag = 2
BarHomogeneous = 1
Presolve = 2
Method = 2 (Method 2 always fails with error code 139 and method 1 takes very long)
MIPGap = 0.1

PreSparsify = 1
Aggregate = 2

Does anyone have any suggestions on how I can solve this? I am happy about all inputs.

Model details
- Framework: Pyomo-based MILP single-objective optimization
- Solver: Gurobi 11.0, academic license
- Variables: 752,720 continuous, 1,046 binary
- Stages: 4 (2020, 2030, 2040, 2050)
- Timesteps: 110 (TSAM-clustered from 8760 hourly data)
- Power unit: GW, Mass unit: tonnes, Currency: kCHF

Last log:

Optimize a model with 1929840 rows, 753766 columns and 4663037 nonzeros
Model fingerprint: 0x8ee96641
Variable types: 752720 continuous, 1046 integer (1046 binary)
Coefficient statistics:
Matrix range [6e-11, 2e+08]
Objective range [1e+00, 1e+00]
Bounds range [5e-03, 2e+08]
RHS range [7e-07, 1e+09]
Warning: Model contains large matrix coefficient range
Presolve removed 765621 rows and 123062 columns (presolve time = 5s) ...
Presolve removed 776898 rows and 124051 columns (presolve time = 10s) ...
Presolve removed 778056 rows and 124075 columns (presolve time = 15s) ...
Presolve removed 778076 rows and 124075 columns (presolve time = 20s) ...
Presolve removed 778076 rows and 124075 columns (presolve time = 25s) ...
Presolve removed 778076 rows and 124075 columns (presolve time = 31s) ...
Presolve removed 778076 rows and 124075 columns (presolve time = 35s) ...
Presolve removed 778076 rows and 124076 columns (presolve time = 40s) ...
Sparsify removed 1534 nonzeros (0%)
Presolve removed 778080 rows and 124063 columns
Presolve time: 42.42s
Presolved: 1151760 rows, 629703 columns, 3577164 nonzeros
Variable types: 629593 continuous, 110 integer (110 binary)
Root relaxation presolve removed 4859 rows and 1423 columns (presolve time = 5s) ...
Root relaxation presolve removed 4859 rows and 1423 columns (presolve time 47 5.81629868e+08 -1.71410883e+09 1.70e+01 4.63e-09 1.26e+03 3746s
resolve time = 15s) ...
Root relaxation presolve removed 9917 rows and 2611 columns (presolve time = 20s) ...
Root relaxation presolve removed 12071 rows and 2946 columns
Root relaxation presolved: 1139689 rows, 626757 columns, 3527133 nonzeros

Root barrier log...

Elapsed ordering time = 7s
Elapsed ordering time = 10s
Elapsed ordering time = 25s
48 5.56695579e+08 -1.51095298e+09 1.55e+01 4.11e-09 1.14e+03 3805s
Ordering time: 135.49s
49 5.12530748e+08 -1.20643385e+09 1.31e+01 2.92e-09 9.47e+02 3869s

Riley Clement

Gurobi Staff

2026年05月13日 00:22

Hi Xinge,

I would not set Aggregate=2, this will be making the problem worse. My guess is that it would be better to also leave PreSparsify at its default. Your other parameter settings look reasonable although ideally would be backed by evidence on smaller models.

You are correct regarding coefficient stats - 19 orders of magnitude is really bad and I am not surprised that Gurobi is struggling to solve this. I would recommend reading the following section in our reference manual for a guide on how to address this: https://docs.gurobi.com/projects/optimizer/en/current/concepts/numericguide/tolerances_scaling.html#advanced-user-scaling

Are you using Python? If you are not sure where the very large, and very tiny coefficients are coming from then you will likely need to do some model analysis, and this is easiest done in Python. I can give you code if so.

We also have an “alpha release” of some scaling functionality in our model-analyzer package:
https://gurobi-modelanalyzer.readthedocs.io/en/latest/scaling.html

Since this is an alpha release, it has not yet been published to PyPI. If you want to try it you will need to install it directly from github:

pip install git+https://github.com/Gurobi/gurobi-modelanalyzer

- Riley

Xinge Gao

2026年05月17日 13:49

Hi Riley,

Thank you for your quick reply and all the information.

I am using Python and would be happy to get the code for finding the small and large coefficients in the matrix range.

Also, I am not fully sure I understand the function of Threads correctly. Is it correct to assume a higher number of Threads (e.g. 16) leads to faster solve times for barrier method? Is there a rule of thumb which number of threads should be used for large-scale energy optimization?

Kindly,

Xinge

2026年05月25日 00:01

You can use the following code

import gurobipy as gp
import numpy as np
import pandas as pd

m = gp.read("/path/to/model")

col_names = {v.index:v.varname for v in m.getVars()}
row_names = {c.index:c.constrname for c in m.getConstrs()}

A = m.getA().tocoo()
df = pd.DataFrame(
    {
        "row": A.row,
        "col": A.col,
        "mag": np.floor(np.log10(abs(A.data))),
    }
)

column_coefficients = (
    df.groupby("col")["mag"]
    .agg(["min", "max"])
    .eval("range = max-min")
    .eval("name = col.map(@col_names)")
    .set_index("name")
)

row_coefficients = (
    df.groupby("row")["mag"]
    .agg(["min", "max"])
    .eval("range = max-min")
    .eval("name = row.map(@row_names)")
    .set_index("name")
)

The “range” column of the column_coefficients and row_coefficients dataframes will be the most useful information. Any value above 6 is not ideal.

You can also use “min” and “max” to find the variables and constraints with small and large coefficients.

Although barrier often does scale well with parallel threads, it often does so up to a point, and then further threads deteriorates the performance. Gurobi may not choose to use all the threads which are available, and will report the number it uses in the “Barrier statistics” section of the log. If it chooses less threads than the number of CPUs on your machine then adding more threads is expected to not help.

Extreme matrix range [6e-11, 2e+08] in multi-stage energy system MILP — model not solving within 24 hours

コメント

お探しのものを見つけられませんでしたか？