Subset constants break associativity: a + factor + b ≠ a + b + factor

[FBumann]

Issue Description

When a constant has a different size on a dimension it shares with a variable/expression, linopy left-joins it — reindexing onto the other operand and dropping whatever does not fit. Because the join drops coordinates, the result depends on the order the terms are written: addition is no longer associative.

Reproducible Example

import pandas as pd
import xarray as xr
import linopy

m = linopy.Model()
a = m.add_variables(coords=[pd.RangeIndex(3, name="time")])  # time 0..2
b = m.add_variables(coords=[pd.RangeIndex(5, name="time")])  # time 0..4
factor = xr.DataArray([10, 20, 30, 40, 50], dims=["time"], coords={"time": range(5)})

print((a + factor + b).const.sel(time=3).item())  # 0.0  — factor dropped
print((a + b + factor).const.sel(time=3).item())  # 40.0 — factor kept

a + factor + b left-joins factor onto a (time 0..2) first, losing time 3–4; a + b + factor joins a + b (time 0..4) first and keeps it.

Expected Behavior

Addition is associative — the result must not depend on term order. Either every order agrees, or the size mismatch raises.

Confirmed on linopy 0.7.0. Part of the arithmetic-convention work (#591).