PyPSA / linopy

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

Re-solving a changed model with "direct" api results in wrong model being solved #266

Closed apfelix closed 1 month ago

apfelix commented 3 months ago

Hi it's me again ^^

I found some strange behaviour using the direct api for highs in resolves. Example: Two simple models. The 1st is max x s.t. x+y <= 1, x>=0, y>=0 The 2nd is max y s.t. x+y <= 1, x>=0.5, y>=0

from linopy import LESS_EQUAL, GREATER_EQUAL, Model

m = Model(chunk=None)

x = m.add_variables(name="x", lower=0.0)
y = m.add_variables(name="y", lower=0.0)

m.add_constraints(x + y, LESS_EQUAL, 1)

m.add_objective(x.to_linexpr(), sense="max")

solver_api=("highs", "lp")
status, condition = m.solve(solver_name=solver_api[0], io_api=solver_api[1])
assert condition == "optimal"
print(m.objective.value)
print(m.solution)

m.add_constraints(x, GREATER_EQUAL, 0.5)
m.add_objective(y.to_linexpr(), sense="max", overwrite=True)

status, condition = m.solve(solver_name=solver_api[0], io_api=solver_api[1])
assert condition == "optimal"
print(m.objective.value)
print(m.solution)

When calling this I get as 2nd solution x = y = 0.5. However, when using the solver_api=("highs", "direct"), I get x = 0.0 and y = 1.0.

I printed the highs models as .lp files after creation and the new constraint x, GREATER_EQUAL, 0.5 did not appear in the 2nd model.

After some digging I found that this is caused by using the @cached_property decorator in the MatrixAccessor class in matrices.py. The model keeps the same MatrixAccessor, stores the values used for building the first model version, and returns them also for the second model version, even though the model has changed. When using @property for all the functions, the correct model is build for me.

I didn't provide a PR for this, as I'm not sure what's the best approach for solving this. Two solutions that I see are that you could either just use regular @property decorators or clear the cache of the corresponding functions after building the model by calling del model.matrices.<cached function name>.

I guess this could affect direct solves with gurobi as well, since similar data structures are used there.

FabianHofmann commented 2 months ago

thanks for raising this @apfelix! it seems to me that deleting the cache after solving makes most sense.