Add support of xarray's diff function to Variable class

[lumbric]

The class linopy.variables.Variable inherits methods from its parent xarray.DataArray. Some of these methods seem to lead to wrong optimization results.

Example:

m = linopy.Model()
var = m.add_variables(name="var", lower=xr.DataArray(np.zeros(5), dims='x', coords={'x': np.arange(5)}))
m.add_constraints(1. * var.diff(dim='x') == 1.)
m.add_objective(var.sum())

Solving using Gurobi gives the solution: var: 0.0 1.0 0.0 0.0 0.0.
The expected solution would be: var: 0.0 1.0 2.0 3.0 4.0

[lumbric]

This does not work either:

m.add_constraints(var.isel(x=slice(None, -1)) - var.isel(x=slice(1, None)) == 1.)

Gurobi fails with "infeasible", the lp file looks like this:

min

obj:

+1.000000 x0
+1.000000 x1
+1.000000 x2
+1.000000 x3
+1.000000 x4

s.t.


c0:
+1.000000 x0
=
1.000000

c1:
+1.000000 x1
-1.000000 x1
=
1.000000

c2:
+1.000000 x2
-1.000000 x2
=
1.000000

c3:
+1.000000 x3
-1.000000 x3
=
1.000000

c4:
-1.000000 x4
=
1.000000


bounds

+0.000000 <= x0 <= +inf
+0.000000 <= x1 <= +inf
+0.000000 <= x2 <= +inf
+0.000000 <= x3 <= +inf
+0.000000 <= x4 <= +inf


binary

end

Using a rule I was able to add the desired constraints:

def diff(m, i):
    if i == 0:
        # workaround, because rule has to return an AnonymousScalarConstraint
        return var[0] - var[0] == 0.
    else:
        return var[i] - var[i-1] == 1.

m.add_constraints(diff, coords=(pd.RangeIndex(5),))

However, I had to add a workaround to omit the constraint for i=0. If there is a better solution, I would be happy to learn! Thanks!

[FabianHofmann]

Hey @lumbric, indeed the diff function should be addressed! You could try out the shift approach instead of substracting (xarray aligns coordinates that's why your second approach was not working). Something like:

m = linopy.Model()
var = m.add_variables(name="var", lower=xr.DataArray(np.zeros(5), dims='x', coords={'x': np.arange(5)}))
lhs = (var - var.shift(x=1)).sel(x=slice(1, None))
m.add_constraints(lhs == 1.)
m.add_objective(var.sum())

[lumbric]

Thanks a lot! This seems to be way faster than the rule-based approach above (2.5s compared to 4.4s for N=100_000 instead of N=5). I don't fully understand why .shift() works and .isel() does not work (see my comment above), but .shift() seems to suite my needs perfectly. Thank you!

(I'll keep the ticket open for the diff() method, feel free to close if you want to track this otherwise or consider it as won't fix.)

[lumbric]

I had a quick look into the implementation, don't really fully grasp all the details yet, basically DataArray methods work if they shift things around (like roll()) or if one of methods overwritten in the Variable class is used to create a LinearExpression, right?
If this is correct, many of the methods inherited from DataArray should not be used then. Maybe some whitelisting mechanism would be a good idea then? Maybe inheritance is not the right pattern here, but something like an adapter? I assume the same holds for the class Constraint and LinearExpression.

I had a very quick look at DataArray methods:

Probably supported correctly by linopy: astype, bfill, ffill, fillna, clip, shift, roll, rolling, where, sum, __getitem__, __repr__, _repr_html_, __neg__, __mul__, __rmul__, __add__, __sub__, __le__, __ge__, __eq__,

Potentially dangerous: diff, T, round, std, __abs__, __and__, __gt__, sortby, transpose, __float__, sel, prod, head, __contains__, __div__, all, any, argmax, argmin, argsort, resample, copy, loc, isel, isin, isnull, groupby, groupby_bins, __floordiv__, __lt__, cumprod, cumsum, differentiate, max, mean, median, min

Harmless, nice to have / necessary: dims, values, size, sizes, shape, name, dtype, __class__, __dir__, __doc__, __format__, attrs, __getattr__, __getattribute__

Harmless, but probably not really needed: to_cdms2, to_dataframe, to_dataset, to_dict, to_index, to_iris, to_masked_array, to_netcdf, to_pandas, to_series

Don't know: to_unstacked_dataset, __enter__, __exit__, drop, drop_sel, drop_vars, dropna, _HANDLED_TYPES, __annotations__, __array__, __array_priority__, __array_ufunc__, __array_wrap__, __bool__, nbytes, plot, rename, chunk, chunks, __complex__, __copy__, __dask_graph__, __dask_keys__, __dask_layers__, __dask_optimize__, __dask_postcompute__, __dask_postpersist__, __dask_scheduler__, __dask_tokenize__, __deepcopy__, __delattr__, __delitem__, __hash__

More methods, didn't look through them yet: __iadd__, __iand__, __ifloordiv__, __imod__, __imul__, __init__, __init_subclass__, __int__, __invert__, __ior__, __ipow__, __isub__, __iter__, __itruediv__, __ixor__, __len__, __matmul__, __mod__, __module__, __ne__, __new__, __or__, __pos__, __pow__, __radd__, __rand__, __reduce__, __reduce_ex__, __rfloordiv__, __rmatmul__, __rmod__, __ror__, __rpow__, __rsub__, __rtruediv__, __rxor__, __setattr__, __setitem__, __sizeof__, __slots__, __str__, __subclasshook__, __truediv__, __weakref__, __xor__, assign_attrs, assign_coords, broadcast_equals, broadcast_like, close, coarsen, combine_first, compute, conj, conjugate, coords, count, data, dot, dt, encoding, equals, expand_dims, from_cdms2, from_dict, from_iris, from_series, get_axis_num, get_index, identical, imag, indexes, integrate, interp, interp_like, interpolate_na, item, load, map_blocks, ndim, notnull, persist, pipe, quantile, rank, real, reduce, reindex, reindex_like, reorder_levels, reset_coords, reset_index, rolling_exp, searchsorted, set_index, squeeze, stack, str, swap_dims, tail, thin, unify_chunks, unstack, var, variable

[FabianHofmann]

Hey @lumbric this is a very good starting point! And I would love to tackle this issue since it would make linopy much safer to use. For Variable and LinearExpression I implemented a decorator which wraps xarray functions, but it is not applied to all relevant functions and does not always makes sense. For LinearExpression I started to disable functions like this but this is incomplete as well.
In the next weeks I won't be able to look into this, but I am always happy to review a PR if you have some spare time (also tackling partially is welcome).

[lumbric]

Okay, will try to see if I find time for a draft

These are possible approaches I can think of at the moment:

• blacklist dangerous methods and attributes as done in LinearExpression for some methods: not the safest way - new xarray attributes are not handled
• whitelist allowed methods using some magic via __getattr__ or so: might involve a lot of magic, not sure if this is really a good idea, but maybe the best compromise
• change the design pattern, do not directly inherit from xarray: would be a major code change, potentially breaking interfaces, I guess)

Do you have any other idea or any opinion on this? I think the most practical solution might be the second one. I can try to make a draft to check if there are any unforeseen downsides of this approach.