PyPSA / linopy

Linear optimization with N-D labeled arrays in Python
https://linopy.readthedocs.io
MIT License
154 stars 42 forks source link

Inconsistent coordinate override in arithmetic expressions #257

Open apfelix opened 3 months ago

apfelix commented 3 months ago

Based on the discussion in #256, I started using <variable>.loc[indices] in arithmetic expressions and observed some inconsistencies (or I still do not understand the behaviour to 100% ^^)

Following setup:

import xarray as xr

import linopy

m = linopy.Model()

# 2 days, hourly (48 hours)
t = [0, 1, 2, 3]
n = ["a", "b"]

level = m.add_variables(coords=[t, n], name="level", dims=["time", "node"])

alpha = xr.DataArray([[1, 2], [3, 4], [5, 6], [7, 8]], coords=level.coords)

other_times = [2, 0, 3, 1]
other_time_level = level.loc[dict(time=other_times)]

Some arithmetic expressions level + other_time_level works as expected, coordinates are overridden based on level

LinearExpression (time: 4, node: 2):
------------------------------------
[0, a]: +1 level[0, a] + 1 level[2, a]
[0, b]: +1 level[0, b] + 1 level[2, b]
[1, a]: +1 level[1, a] + 1 level[0, a]
[1, b]: +1 level[1, b] + 1 level[0, b]
[2, a]: +1 level[2, a] + 1 level[3, a]
[2, b]: +1 level[2, b] + 1 level[3, b]
[3, a]: +1 level[3, a] + 1 level[1, a]
[3, b]: +1 level[3, b] + 1 level[1, b]

When using level + alpha * other_time_level, I would have expected to get

LinearExpression (time: 4, node: 2):
------------------------------------
[0, a]: +1 level[0, a] + 1 level[2, a]
[0, b]: +1 level[0, b] + 2 level[2, b]
[1, a]: +1 level[1, a] + 3 level[0, a]
[1, b]: +1 level[1, b] + 4 level[0, b]
[2, a]: +1 level[2, a] + 5 level[3, a]
[2, b]: +1 level[2, b] + 6 level[3, b]
[3, a]: +1 level[3, a] + 7 level[1, a]
[3, b]: +1 level[3, b] + 8 level[1, b]

since I though that alpha would be term to define the coordinates in alpha * other_time_level

However, what I actually get is the following (coordinates of other_time_level are sorted according to the original time index again)

LinearExpression (time: 4, node: 2):
------------------------------------
[0, a]: +1 level[0, a] + 1 level[0, a]
[0, b]: +1 level[0, b] + 2 level[0, b]
[1, a]: +1 level[1, a] + 3 level[1, a]
[1, b]: +1 level[1, b] + 4 level[1, b]
[2, a]: +1 level[2, a] + 5 level[2, a]
[2, b]: +1 level[2, b] + 6 level[2, b]
[3, a]: +1 level[3, a] + 7 level[3, a]
[3, b]: +1 level[3, b] + 8 level[3, b]

Interestingly, the result changes when transforming other_time_level into a LinearExpression first, i.e., call level + alpha * other_time_level.to_linexpr()

LinearExpression (time: 4, node: 2):
------------------------------------
[0, a]: +1 level[0, a] + 5 level[2, a]
[0, b]: +1 level[0, b] + 6 level[2, b]
[1, a]: +1 level[1, a] + 1 level[0, a]
[1, b]: +1 level[1, b] + 2 level[0, b]
[2, a]: +1 level[2, a] + 7 level[3, a]
[2, b]: +1 level[2, b] + 8 level[3, b]
[3, a]: +1 level[3, a] + 3 level[1, a]
[3, b]: +1 level[3, b] + 4 level[1, b]

Finally, I tried to simply add alpha instead of multiplying it, i.e., level + alpha + other_time_level and this gives the expected result in terms of sorting the single terms into the expected rows:

LinearExpression (time: 4, node: 2):
------------------------------------
[0, a]: +1 level[0, a] + 1 level[2, a] + 1
[0, b]: +1 level[0, b] + 1 level[2, b] + 2
[1, a]: +1 level[1, a] + 1 level[0, a] + 3
[1, b]: +1 level[1, b] + 1 level[0, b] + 4
[2, a]: +1 level[2, a] + 1 level[3, a] + 5
[2, b]: +1 level[2, b] + 1 level[3, b] + 6
[3, a]: +1 level[3, a] + 1 level[1, a] + 7
[3, b]: +1 level[3, b] + 1 level[1, b] + 8

So my question would be: How to achieve the intended result and is all of the behaviour as it should be?

FabianHofmann commented 3 months ago

very good that you look into this @apfelix, there is some inconsistency indeed. However, I would not fully have the same expectation on the operations.

  1. In my point of view overriding should only happen when adding and subtracting terms. works as intended. See level + other_time_level

  2. multiplication should be aligned by coordinates, no overriding. Thus, when multiplying a variable with a coefficient, the first object in the multiplication should determine the coordinate alignment of the result . Meaning:

    alpha * other_time_level -> result has the same coordinates as alpha, other_time_level is sorted accordingly. works as intended

    other_time_level * alpha -> result has the same coordinates as other_time_level, alpha is sorted accordingly. :negative_squared_cross_mark: does not work.

  3. If one wants to ignore the coefficient coords in the multiplication, i.e. not aligning the coefficient with the variable, then one should use

    other_time_level * alpha.values

what do you think?

apfelix commented 3 months ago

I though a little about this, and I think your expectation does make sense. However, I'm curious on what should happen if the multiplication does involve two variables. Currently level * other_time_level also works as coefficient multiplication, i.e., aligning the values by coordinate again

QuadraticExpression (time: 4, node: 2):
---------------------------------------
[0, a]: +1 level[0, a] level[0, a]
[0, b]: +1 level[0, b] level[0, b]
[1, a]: +1 level[1, a] level[1, a]
[1, b]: +1 level[1, b] level[1, b]
[2, a]: +1 level[2, a] level[2, a]
[2, b]: +1 level[2, b] level[2, b]
[3, a]: +1 level[3, a] level[3, a]
[3, b]: +1 level[3, b] level[3, b]

But then how would you achieve

QuadraticExpression (time: 4, node: 2):
---------------------------------------
[0, a]: +1 level[0, a] level[2, a]
[0, b]: +1 level[0, b] level[2, b]
[1, a]: +1 level[1, a] level[0, a]
[1, b]: +1 level[1, b] level[0, b]
[2, a]: +1 level[2, a] level[3, a]
[2, b]: +1 level[2, b] level[3, b]
[3, a]: +1 level[3, a] level[1, a]
[3, b]: +1 level[3, b] level[1, b]

?

If the multiplication of linopy objects should also ignore alignment, my suggestion for a consistent set of alignment rules would be the following (not sure how hard it would be to implement this):

Order of the resulting coordinates should always be defined by the first term that actually has coordinates, so numpy arrays would be ignored even if they are first.

FabianHofmann commented 3 months ago

Interesting, let me think about it. I would like to have the API change as small as possible, as changes can lead to very much unexpected behaviour.