Fix timeslice-related bugs

[tsmbland]

This fixes a few subtle bugs which were very difficult to spot, but messing up some of the results. It's all to do with timeslices and broadcasting...

---

There are many types of variable in MUSE (supply, demand, capacity, costs, technodata variables etc.), some of which may differ between timeslices (e.g. supply) and some of which are constant for the year (e.g. capacity). MUSE uses xarray dataarrays/datasets to represent these, which will either have a timeslice dimension (multi-index) or not depending on whether the data is timesliced.

In a timesliced dataset, values for each timeslice correspond to the time-period of the timeslice. For example, if you had a timesliced supply dataset, the value for each timeslice would correspond to the total supply over the length of the timeslice. If you had a supply dataset without a timeslice dimension, the value would correspond to supply across the year (i.e. the sum across timeslices).

You can use the convert_timeslice function to extend a non-timesliced dataset over the timeslice dimension. There are a couple of ways you can do this, depending on the variable (although see #516). For example, if you had supply data for the year, you could split this up so that each timeslice has a proportional supply value according to its length, so that the sum across timeslices is equal to the original yearly value. Conversely, if you had prices data, you'd want to copy/broadcast your yearly price across all timeslices rather than splitting it up.

That's pretty much the background, which is all straightforward enough. However, in practice, there are lots of subtleties, and it's very easy to make mistakes when dealing with a mix of timesliced and non-timesliced objects. This PR addresses three such mistakes, which seem to be responsible for the issue described in #512

---

The first is to do with the utilization factor. The UF can either be specified by the user for the year as a whole or separately for each timeslice, and the internal dataarray will either have a timeslice dimension or not depending on how the UF is specified in the input data. Unfortunately, this means you have to be very careful with some of the maths. Let's consider a simple example with four timeslices, an output rate of 1/year and a utilization factor of 1 in all timeslices. We want to calculate the maximum output in each timeslice, which is simply the output rate for the timeslice multiplied by the utilization factor for the timeslice. We can use convert_timeslice to split the yearly output rate over the four timeslices, but then how should we deal with the UF?

This is how it's currently coded:
max_production = convert_timeslice(fixed_outputs * utilization_factor)

But this means that you end up with a different answer depending on whether the utilization factor is specified for the year as a whole or for each timeslice (even if each timeslice has the same UF):

With a single utilization factor:

convert_timeslice(1 * 1)
-> convert_timeslice(1)
-> [0.25, 0.25, 0.25, 0.25]

With timeslice-level utilization factor:

convert_timeslice(1 * [1, 1, 1, 1])
-> convert_timeslice([1, 1, 1, 1])
-> [1, 1, 1, 1]

However, if we change this to the following it becomes more consistent
max_production = convert_timeslice(fixed_outputs) * utilization_factor

With a single utilization factor:

convert_timeslice(1) * 1
-> [0.25, 0.25, 0.25, 0.25] * 1
-> [0.25, 0.25, 0.25, 0.25]

With timeslice-level utilization factor:

convert_timeslice(1) * [1, 1, 1, 1]
-> [0.25, 0.25, 0.25, 0.25] *  [1, 1, 1, 1]
-> [0.25, 0.25, 0.25, 0.25]

Very subtle, and would never be spotted without digging into the code in debug mode, but this can make a big difference to the results.

---

The second is to do with the maximum_production function, which would sometimes return timesliced data and sometimes aggregated yearly data depending on the context in which it was called. Again, you have to be very careful with the maths when using this data.

Let's say you're subtracting production from demand to calculate excess demand. If your demand data is timesliced and your production data isn't, then the subtraction operation will broadcast your production data across the timeslices before subtracting from demand. This means you end up subtracting the full yearly production from each timesliced demand figure, which obviously isn't appropriate. You could use convert_timeslice on the production data beforehand, but I think the better approach is the ensure that max_production always returns timesliced data, which I've done here. The net result of this was to fix an issue with the _inner_split function like what I've just described (which seems to manifest only when not using retrofit agents). Again, near impossible to spot without some serious debugging.

---

The third issue is to do with the ordering of timeslices. When the technodata_timeslices csv is loaded (contains timeslice-level UF and MSF data), the timeslices are sorted into alphabetical order in the resulting dataset, whereas all other timesliced objects match the order specified in the settings file. This isn't a problem when you're performing operations on the xarray object like above, as xarray will use timeslice names to match datasets. But this may be a problem when preparing numpy arrays for the solver, as the ordering in these arrays will match the ordering of the xarray object it's derived from.

I've adjusted the read_technodata_timeslices function to sort the object when it's read, which should clear up any inconsistencies.

---

Overall, the issue described in #512 appears to be fixed.

However, given how complicated and difficult to spot these bugs were, I'm not completely convinced that I've fixed everything. I think we need a more radical approach to guarantee that these kinds of errors can never exist.

A good start would be to decide, for each variable, whether it should exist in timesliced form or yearly form, and make sure this is consistent everywhere in the code regardless of the input data. For example:

• supply/demand: should always be timesliced
• capacity: should always be yearly
• commodity prices: should always be timesliced
• technology variables apart from UF/MSF: should always be yearly
• UF/MSF: should always be timesliced

I'm not really sure how best to do this, apart from scattering a load of assert statements throughout the code which could be a bit messy.

The codebase is also far looser than it needs to be in terms of potentially allowing different objects to have different timeslicing schemes. I don't think there will ever be a need for this (except maybe the legacy sectors), so I think we should make use of the global TIMESLICE object (which is set according to what's specified in the settings file), and make sure that all timesliced objects match this scheme. This is the basis of #519, although that's still a work in progress.

I'd be a lot happier if automatic broadcasting across the timeslice dimension was banned, because this can lead to some very hard to spot bugs. But I'm not really sure how we can do that.

---

The results have changed for all models, but the changes are mostly very small. I have also modified the min/max timeslice tutorial to make it clearer what the expected results are supposed to be

You can see the tutorial notebooks here

Closes #512

[tsmbland]

@dalonsoa Any ideas how to make the code more robust so that this sort of thing can never happen again? It's really hard to spot these kinds of bugs just by looking at the code, and obviously the tests good enough to catch them...

[dalonsoa]

These changes all look sensible to me, therefore approving. But, as you said, there might be plenty of other places where the use of timeslices and broadcasting is inconsistent - or plainly wrong - and very difficult to catch.

There was an attempt to add a arithmetic_broadcast=False global flag to prevent automatic broasdcasting. Unfortunately, there were some issues related to dask (odd issues) and the change was reverted.

I'm not sure how good practice this would be or if it would work, but you could implement the changes in the above PR - some of them - as a patch. Something along the lines of

# In muse/__main__.py

def patched_broadcast_compat_data(self, other):
    from xarray.core.variable import Variable

    if (isinstance(other, Variable) and self.dims != other.dims) or (
            is_duck_array(other) and self.ndim != other.ndim
        ):
            raise ValueError(
                "Broadcasting is necessary but automatic broadcasting is disabled globaly."
            )

    if all(hasattr(other, attr) for attr in ["dims", "data", "shape", "encoding"]):
        # `other` satisfies the necessary Variable API for broadcast_variables
        new_self, new_other = _broadcast_compat_variables(self, other)
        self_data = new_self.data
        other_data = new_other.data
        dims = new_self.dims
    else:
        # rely on numpy broadcasting rules
        self_data = self.data
        other_data = other
        dims = self.dims
    return self_data, other_data, dims

...

if "__main__" == __name__:
    from unittest.mock import patch
    with patch("xarray.core.variable._broadcast_compat_data", patched_broadcast_compat_data):
        run()

Assuming it works, it won't affect the tests, but it should pick any attempt to do automatic broadcasting when actually running a model, which is what we want, in the end.

As I said, I've no idea if this would work, but it should be possible to do something along these lines.

Also, obviously, include the explanation you give in the PR in a developers documentation section.

[dalonsoa]

I think this is what confuses me the most about timeslices, specially now that you have shed some light into it: when do you use the timeslice of one particular array, like here, and when the global TIMESLICE? If they are the same I'd use always the global one, for clarity. And if they are not... why they are not and how can we know?

[tsmbland]

I'm pretty sure these are all equivalent to TIMESLICE (or at least they should be), so there's no reason not to use the global. This is part of what I'm doing in #519, although that PR has become a bit too big so I might try to break this specific change into its own PR

[tsmbland]

Actually it's more complicated than this. TIMESLICE is the global timeslicing scheme from the timeslices section of the settings file. However, MUSE does allow you to have different timeslicing schemes for different sectors (see here). I imagine this is so you can have less granularity in some sectors. For example, in the oil sector, you may only care about meeting demands at the seasonal level, so your timeslices might be "winter", "summer" etc. rather than "winter.weekend.morming", "winter.weekday.night" etc. In this case, maybe the timeslicing of the array isn't going to match up with TIMESLICE. Not sure whether it actually works like that though...

In that case, #519 is a waste of time as it gets rid of this functionality, although I do still want to tidy up how timeslices are dealt with as it's a complete mess at the moment.

[tsmbland]

Thanks! This is really useful. I'm going to give this a try and see how I get on