Memory error for transport model using scipy solver

[tsmbland]

I tried using the scipy solver with @sharwaridixit's transport model, and got this rather alarming error:

numpy.core._exceptions._ArrayMemoryError: Unable to allocate 40.2 GiB for an array with shape (16, 16, 7, 7, 41, 41, 16, 16) and data type float64

It's specifically being thrown by this line in the calculation of the constraint matrix. I'm guessing there's a bug somewhere that's accidentally duplicating all the dimensions of an array, but I'd need to investigate further

[HarmonicReflux]

Hi @tsmbland,

Can you let me know what command you invoke in the terminal that produces the error, i.e. is it something like
muse --model -<your_model_goes_here>

I take from your description, the culprit is line 853 in src/muse/constraints.py and reads
reduced = reduce(xr.DataArray.__mul__, diagonal_submats)

[tsmbland]

Hi @HarmonicReflux. To reproduce the error you will need to use the model found here. However, in order to fix the bug described in that issue and get it working with the latest version of MUSE, you will have to replace the settings file with the one attached here:
settings.zip
You can then run the model with muse settings.toml.

The specific issue I'm describing here only happens when you set lpsolver to "scipy" (currently it's using "adhoc").

Let me know if you have any problems. I imagine this one won't be super straightforward to fix.

[HarmonicReflux]

• How can one access the source code for the LP solver "adhoc" and where is this solver documented?

• I think the solver "adhoc" acts fundamentally different from scipy's version of solving a linear program, although there
  are several solvers available for LP's and settings.toml does not specify which one it uses. In particular there are "highs", "highs-ds", "interior-point" (legacy), "revised simplex" (legacy), and "simplex" (legacy) available.

In the below, I print a detailed error log:

My take is the scipy solver simply runs out of memory for the (16, 16, 7, 7, 41, 41, 16, 16) data structure, which has about 5.4 billion entries:.

When we decrease the problem to solve, for example, within settings.toml, shrink

time_framework = [2010, 2015, 2020, 2025, 2030, 2035, 2040, 2045, 2050]
foresight = 5  

to

time_framework = [2010, 2020, 2030, 2040, 2050]
foresight = 10  # Multiple of the smallest difference

we get a numpy.core._exceptions._ArrayMemoryError: Unable to allocate 30.8 GiB for an array with shape (14, 14, 7, 7, 41, 41, 16, 16) and data type float64. This means are now down to 31 GB from the original 40GB (although, we then get a ValueError saying `File "C:\Users\bjs13\python_venvs\muse_os\lib\site-packages\xarray\core\variable.py", line 3001, in calculate_dimensions
raise ValueError(
ValueError: dimension 'year' already exists as a scalar variable
should we chose the "adhoc" solver on the smaller problem.

There are also commercial solvers for LP's such as "Gurobi" and "CPLex" which have free academic licenses, and may be better suited to solve the problem?

[ahawkes]

On the memory error, what is this massive thing with 3.5 billion entries. Seems ludicrously big for the problem at hand?

…

________________________________ From: HarmonicReflux ***@***.***> Sent: Monday, July 15, 2024 4:46:56 PM To: EnergySystemsModellingLab/MUSE_OS ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [EnergySystemsModellingLab/MUSE_OS] Memory error for transport model using scipy solver (Issue #389) This email from ***@***.*** originates from outside Imperial. Do not click on links and attachments unless you recognise the sender. If you trust the sender, add them to your safe senders list<https://spam.ic.ac.uk/SpamConsole/Senders.aspx> to disable email stamping for this address. * How can one access the source code for the LP solver "adhoc" and where is this solver documented? * I think the solver "adhoc" acts fundamentally different from scipy's version of solving a linear program, although there are several solvers available <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html> for LP's and settings.toml does not specify which one it uses. In particular there are "highs", "highs-ds", "interior-point" (legacy), "revised simplex" (legacy), and "simplex" (legacy) available. In the below, I print a detailed error log: lp_solver_scipy.jpg (view on web)<https://github.com/user-attachments/assets/31fb0d9c-c66a-4e0b-a8ed-0f6d71d8e420> My take is the scipy solver simply runs out of ram for the (16, 16, 7, 7, 41, 41, 16, 16) data structure, which has about 5.4 billion entries:. When we decrease the problem to solve, for example, within settings.toml, shrink time_framework = [2010, 2015, 2020, 2025, 2030, 2035, 2040, 2045, 2050] foresight = 5 to time_framework = [2010, 2020, 2030, 2040, 2050] foresight = 10 # Multiple of the smallest difference we get a numpy.core._exceptions._ArrayMemoryError: Unable to allocate 30.8 GiB for an array with shape (14, 14, 7, 7, 41, 41, 16, 16) and data type float64. This means are now down to 31 GB from the original 42GB (although, we then get a ValueError saying `File "C:\Users\bjs13\python_venvs\muse_os\lib\site-packages\xarray\core\variable.py", line 3001, in calculate_dimensions raise ValueError( ValueError: dimension 'year' already exists as a scalar variable should we chose the "adhoc" solver on the smaller problem. There are also commercial solvers for LP's such as "Gurobi" and "CPLex" which have free academic licenses, and may be better suited to solve the problem? — Reply to this email directly, view it on GitHub<#389 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AC37JLLGEGKK4WLAVYEIG53ZMPVHBAVCNFSM6AAAAABKHCCZWGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMRYHAZDIMBTHA>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[HarmonicReflux]

Looks very fishy indeed.

To isolate the issue and compare the behaviour of the optimiser for the output of
muse --model default, I stored them in a their separate folder on my local copy of the repository.

Say, I work from the project root, and can run the adhoc optimiser fine with
muse my_models/Transporttrialmodel/settings.toml

If I now use your settings.toml file, and apply it to the Results of the default model like so:
muse my_models/Defaultmodel/settings.toml

I get the following AssertionError (see also screenshot attached):
AssertionError: Projections file does not exist (/home/bjs/Downloads/vcs/MUSE_OS/my_models/Defaultmodel/input/Projections.csv)

Is there a quick way to adapt the settings.toml file as well as the required input files to make it compliant to the structure of the Results file for the default model?

[HarmonicReflux]

Another idea I had is to experimentally run the optimizer of scipy for all algorithms supported off-the-shelf, i.e.
"highs", which is the default, "highs-ds", "interior-point" (legacy), "revised simplex" (legacy), and "simplex" (legacy)

How can we manipulate the settings.toml file to cater for these different flavours, say, by passing a respective keyword argument onto the solver?

[ahawkes]

I think anything that says “legacy” here relates to an older version of the model (StarMuse). So I guess only highs and highs-ds are relevant. Not sure if the difference between these.

…

________________________________ From: HarmonicReflux ***@***.***> Sent: Thursday, July 18, 2024 9:18:09 AM To: EnergySystemsModellingLab/MUSE_OS ***@***.***> Cc: Hawkes, Adam D ***@***.***>; Comment ***@***.***> Subject: Re: [EnergySystemsModellingLab/MUSE_OS] Memory error for transport model using scipy solver (Issue #389) This email from ***@***.*** originates from outside Imperial. Do not click on links and attachments unless you recognise the sender. If you trust the sender, add them to your safe senders list<https://spam.ic.ac.uk/SpamConsole/Senders.aspx> to disable email stamping for this address. Another idea I had is to experimentally run the optimizer of scipy for all algorithms supported off-the-shelf, i.e. "highs", which is the default, "highs-ds", "interior-point" (legacy), "revised simplex" (legacy), and "simplex" (legacy) How can we manipulate the settings.toml file to cater for these different flavours, say, by passing a respective keyword argument onto the solver? — Reply to this email directly, view it on GitHub<#389 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AC37JLKADIOVJPFUJTDCIXDZM524DAVCNFSM6AAAAABKHCCZWGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMZVHEYDQNBYGU>. You are receiving this because you commented.Message ID: ***@***.***>

[tsmbland]

I've had a look into this a bit more, and it's not a bug per-se, rather a combinatorial explosion due to the large number of timeslices, technologies and end-use commodities (particularly in the road_vehicles sector).

The particular object facing the memory issue is the constraint matrix for the max_production constraint. For this constraint, A is a square matrix with side length (n_assets * n_technologies * n_commodities * n_timeslices), which for this model is 73472 x 73472 (i.e. very big!)

There are a couple of ways we can address this, either:

Changing the model:

• The model has 16 timeslices, although none of the parameters/inputs are defined on a per-timeslice basis, so there's no need to have this many timeslices (or even more than one timeslice)

Changing MUSE:

• The max_production constraint (and probably other constraints) could be compressed over the timeslice dimension if there are no differences between timeslices.
• Also, this constraint is being defined on a per-commodity basis, but there's no reason to do this as max_production is defined in terms of technology flow not per-commodity (indeed the ratio between the output commodities of a technology should be fixed by definition)

This latter point brings me to a bigger issue. The output of the linear optimisation is a decision vector containing capacity and production data. The production data is defined for every commodity, but there's nothing in place to constrain the ratio of commodity outputs for a given technology. Surely, instead of calculating production of each commodity, we should be calculating the flow of each technology, and then from this we can calculate commodity outputs based on fixed output ratios in a separate step.

I also think this is why we cannot currently use the production data directly from the solver, as described in #423