diff --git a/doc/index.rst b/doc/index.rst index fd7f9ed8..3d37e8d2 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -116,6 +116,7 @@ This package is published under MIT license. piecewise-linear-constraints piecewise-linear-constraints-tutorial manipulating-models + iterative-resolving testing-framework transport-tutorial infeasible-model diff --git a/doc/iterative-resolving.nblink b/doc/iterative-resolving.nblink new file mode 100644 index 00000000..1492aeb9 --- /dev/null +++ b/doc/iterative-resolving.nblink @@ -0,0 +1,3 @@ +{ + "path": "../examples/iterative-resolving.ipynb" +} diff --git a/doc/release_notes.rst b/doc/release_notes.rst index b4a92e64..035209e5 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -4,6 +4,8 @@ Release Notes Upcoming Version ---------------- +* Add ``Solver.resolve()`` method for re-solving a modified native solver model without rebuilding the linopy model. Implemented for HiGHS, Gurobi, and Mosek. +* Add ``Model.apply_result()`` to map a solver ``Result`` back onto model variables and constraints, enabling iterative solve workflows. * Harmonize coordinate alignment for operations with subset/superset objects: - Multiplication and division fill missing coords with 0 (variable doesn't participate) - Addition and subtraction of constants fill missing coords with 0 (identity element) and pin result to LHS coords diff --git a/examples/iterative-resolving.ipynb b/examples/iterative-resolving.ipynb new file mode 100644 index 00000000..6ebdb8e1 --- /dev/null +++ b/examples/iterative-resolving.ipynb @@ -0,0 +1,302 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d21785bf", + "metadata": {}, + "source": [ + "# Iterative Re-solving\n", + "\n", + "In many optimization workflows you need to solve a model, tweak parameters on the native solver object, and re-solve — without rebuilding the entire linopy model. This is common in:\n", + "\n", + "- Sensitivity analysis (varying objective coefficients or bounds)\n", + "- Decomposition algorithms (Benders, column generation)\n", + "- Rolling-horizon / receding-horizon schemes\n", + "\n", + "Linopy supports this via two methods:\n", + "\n", + "- **`Solver.resolve()`** — re-solves an existing native solver model and returns a `Result`\n", + "- **`Model.apply_result()`** — maps a `Result` back onto the linopy model's variables and constraints\n", + "\n", + "Currently `resolve()` is implemented for the **HiGHS**, **Gurobi**, and **Mosek** solvers." + ] + }, + { + "cell_type": "markdown", + "id": "f7f604f5", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "We start with a simple two-variable LP." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "122ee88e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-16T08:48:03.207441Z", + "iopub.status.busy": "2026-03-16T08:48:03.207185Z", + "iopub.status.idle": "2026-03-16T08:48:03.649145Z", + "shell.execute_reply": "2026-03-16T08:48:03.648818Z" + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from linopy import Model\n", + "from linopy.solvers import Highs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c05110c", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-16T08:48:03.650662Z", + "iopub.status.busy": "2026-03-16T08:48:03.650440Z", + "iopub.status.idle": "2026-03-16T08:48:03.744258Z", + "shell.execute_reply": "2026-03-16T08:48:03.743824Z" + } + }, + "outputs": [], + "source": [ + "m = Model()\n", + "x = m.add_variables(lower=0, upper=10, name=\"x\")\n", + "y = m.add_variables(lower=0, upper=10, name=\"y\")\n", + "m.add_constraints(x + y >= 8, name=\"demand\")\n", + "m.add_constraints(x + y <= 15, name=\"capacity\")\n", + "m.objective = 2 * x + 3 * y\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "574be9bd", + "metadata": {}, + "source": [ + "## Initial Solve\n", + "\n", + "Solve the model as usual. After solving, linopy stores the native solver object in `m.solver_model`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "548b7a4d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-16T08:48:03.745807Z", + "iopub.status.busy": "2026-03-16T08:48:03.745531Z", + "iopub.status.idle": "2026-03-16T08:48:03.787345Z", + "shell.execute_reply": "2026-03-16T08:48:03.786876Z" + } + }, + "outputs": [], + "source": [ + "m.solve(solver_name=\"highs\")\n", + "print(f\"Objective: {m.objective.value}\")\n", + "print(f\"x = {m.solution['x'].values.item():.2f}\")\n", + "print(f\"y = {m.solution['y'].values.item():.2f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "214cbe09", + "metadata": {}, + "source": [ + "## Modify the Native Solver Model\n", + "\n", + "Access the native HiGHS object and change the objective coefficients. Here we make `x` much more expensive, so the solver should prefer `y`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b662db35", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-16T08:48:03.789208Z", + "iopub.status.busy": "2026-03-16T08:48:03.788835Z", + "iopub.status.idle": "2026-03-16T08:48:03.801031Z", + "shell.execute_reply": "2026-03-16T08:48:03.800672Z" + } + }, + "outputs": [], + "source": [ + "h = m.solver_model\n", + "n_cols = h.getNumCol()\n", + "new_costs = np.array([10.0, 1.0], dtype=float)\n", + "h.changeColsCost(n_cols, np.arange(n_cols, dtype=np.int32), new_costs)" + ] + }, + { + "cell_type": "markdown", + "id": "9e55e60b", + "metadata": {}, + "source": [ + "## Re-solve and Apply\n", + "\n", + "Create a solver instance, call `resolve()`, then apply the result back to the linopy model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "949ce62e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-16T08:48:03.802549Z", + "iopub.status.busy": "2026-03-16T08:48:03.802443Z", + "iopub.status.idle": "2026-03-16T08:48:03.818907Z", + "shell.execute_reply": "2026-03-16T08:48:03.818509Z" + } + }, + "outputs": [], + "source": [ + "solver = Highs()\n", + "result = solver.resolve(h, sense=m.sense)\n", + "m.apply_result(result, solver_name=\"highs\")\n", + "\n", + "print(f\"Objective: {m.objective.value}\")\n", + "print(f\"x = {m.solution['x'].values.item():.2f}\")\n", + "print(f\"y = {m.solution['y'].values.item():.2f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "5b2eab9a", + "metadata": {}, + "source": [ + "Since `x` is now 10× more expensive than `y`, the solver allocates as much as possible to `y`." + ] + }, + { + "cell_type": "markdown", + "id": "7fd50c4d", + "metadata": {}, + "source": [ + "## Iterating Over Parameter Sweeps\n", + "\n", + "The `resolve` / `apply_result` pattern is efficient for parameter sweeps because it avoids rebuilding the model from scratch each time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb3ac524", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-16T08:48:03.820485Z", + "iopub.status.busy": "2026-03-16T08:48:03.820382Z", + "iopub.status.idle": "2026-03-16T08:48:03.862574Z", + "shell.execute_reply": "2026-03-16T08:48:03.862143Z" + } + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "cost_x_values = [1.0, 2.0, 5.0, 10.0, 20.0]\n", + "results = []\n", + "\n", + "for cost_x in cost_x_values:\n", + " h.changeColsCost(\n", + " n_cols, np.arange(n_cols, dtype=np.int32), np.array([cost_x, 3.0], dtype=float)\n", + " )\n", + " result = solver.resolve(h, sense=m.sense)\n", + " m.apply_result(result, solver_name=\"highs\")\n", + " results.append(\n", + " {\n", + " \"cost_x\": cost_x,\n", + " \"objective\": m.objective.value,\n", + " \"x\": m.solution[\"x\"].values.item(),\n", + " \"y\": m.solution[\"y\"].values.item(),\n", + " }\n", + " )\n", + "\n", + "pd.DataFrame(results)" + ] + }, + { + "cell_type": "markdown", + "id": "6f84f855", + "metadata": {}, + "source": [ + "## Handling Infeasible Results\n", + "\n", + "When a re-solve yields an infeasible or otherwise non-optimal result, `apply_result` returns the status without setting variable solutions. You can also pass `solution=None` in a `Result` object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acc2e8bb", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-16T08:48:03.864315Z", + "iopub.status.busy": "2026-03-16T08:48:03.864083Z", + "iopub.status.idle": "2026-03-16T08:48:03.879731Z", + "shell.execute_reply": "2026-03-16T08:48:03.879245Z" + } + }, + "outputs": [], + "source": [ + "from linopy.constants import Result, Status, TerminationCondition\n", + "\n", + "status = Status.from_termination_condition(TerminationCondition.infeasible)\n", + "result = Result(status=status, solution=None)\n", + "\n", + "s, tc = m.apply_result(result)\n", + "print(f\"Status: {s}, Termination: {tc}\")" + ] + }, + { + "cell_type": "markdown", + "id": "898f4824", + "metadata": {}, + "source": [ + "## API Reference\n", + "\n", + "**`Solver.resolve(solver_model, sense=\"min\")`**\n", + "\n", + "Re-solves a native solver object and returns a `Result`. Implemented for `Highs`, `Gurobi`, and `Mosek`. Raises `NotImplementedError` for other solvers.\n", + "\n", + "**`Model.apply_result(result, solver_name=None)`**\n", + "\n", + "Applies a `Result` to the model:\n", + "1. Resets existing solutions\n", + "2. Sets status and termination condition\n", + "3. Maps primal values back to variables\n", + "4. Maps dual values back to constraints (if available)\n", + "\n", + "Returns a `(status, termination_condition)` tuple." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/linopy/model.py b/linopy/model.py index 54334411..ae32ab5b 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -44,6 +44,7 @@ SOS_TYPE_ATTR, TERM_DIM, ModelStatus, + Result, TerminationCondition, ) from linopy.constraints import AnonymousScalarConstraint, Constraint, Constraints @@ -1553,50 +1554,77 @@ def solve( os.remove(fn) try: - result.info() - - self.objective._value = result.solution.objective - self.status = result.status.status.value - self.termination_condition = result.status.termination_condition.value - self.solver_model = result.solver_model - self.solver_name = solver_name - - if not result.status.is_ok: - return ( - result.status.status.value, - result.status.termination_condition.value, - ) + return self.apply_result(result, solver_name) + finally: + if sos_reform_result is not None: + undo_sos_reformulation(self, sos_reform_result) - # map solution and dual to original shape which includes missing values - sol = result.solution.primal.copy() - sol = set_int_index(sol) - sol.loc[-1] = nan + def apply_result( + self, + result: Result, + solver_name: str | None = None, + ) -> tuple[str, str]: + """ + Apply a solver Result to the model, mapping primal/dual values back + to variables and constraints. - for name, var in self.variables.items(): - idx = np.ravel(var.labels) + Useful for iterative re-solve workflows where the solver model is + modified and re-solved outside of ``Model.solve()``. + + Parameters + ---------- + result : Result + Result object from a solver's ``resolve()`` or ``solve_problem_*`` method. + solver_name : str, optional + Name of the solver used. + + Returns + ------- + tuple[str, str] + Status and termination condition. + """ + self.reset_solution() + + self.status = result.status.status.value + self.termination_condition = result.status.termination_condition.value + self.solver_model = result.solver_model + self.solver_name = solver_name or "unknown" + + if result.solution is None: + return (self.status, self.termination_condition) + + result.info() + self.objective._value = result.solution.objective + + if not result.status.is_ok: + return (self.status, self.termination_condition) + + sol = result.solution.primal.copy() + sol = set_int_index(sol) + sol.loc[-1] = nan + + for name, var in self.variables.items(): + idx = np.ravel(var.labels) + try: + vals = sol[idx].values.reshape(var.labels.shape) + except KeyError: + vals = sol.reindex(idx).values.reshape(var.labels.shape) + var.solution = xr.DataArray(vals, var.coords) + + if not result.solution.dual.empty: + dual = result.solution.dual.copy() + dual = set_int_index(dual) + dual.loc[-1] = nan + + for name, con in self.constraints.items(): + idx = np.ravel(con.labels) try: - vals = sol[idx].values.reshape(var.labels.shape) + vals = dual[idx].values.reshape(con.labels.shape) except KeyError: - vals = sol.reindex(idx).values.reshape(var.labels.shape) - var.solution = xr.DataArray(vals, var.coords) - - if not result.solution.dual.empty: - dual = result.solution.dual.copy() - dual = set_int_index(dual) - dual.loc[-1] = nan - - for name, con in self.constraints.items(): - idx = np.ravel(con.labels) - try: - vals = dual[idx].values.reshape(con.labels.shape) - except KeyError: - vals = dual.reindex(idx).values.reshape(con.labels.shape) - con.dual = xr.DataArray(vals, con.labels.coords) - - return result.status.status.value, result.status.termination_condition.value - finally: - if sos_reform_result is not None: - undo_sos_reformulation(self, sos_reform_result) + vals = dual.reindex(idx).values.reshape(con.labels.shape) + con.dual = xr.DataArray(vals, con.labels.coords) + + return self.status, self.termination_condition def _mock_solve( self, diff --git a/linopy/solvers.py b/linopy/solvers.py index 10731547..893f8478 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -19,7 +19,7 @@ from collections import namedtuple from collections.abc import Callable, Generator from pathlib import Path -from typing import TYPE_CHECKING, Any, Generic, TypeVar +from typing import TYPE_CHECKING, Any, Generic, Literal, TypeVar import numpy as np import pandas as pd @@ -419,6 +419,31 @@ def solve_problem( msg = "No problem file or model specified." raise ValueError(msg) + def resolve( + self, solver_model: Any, sense: Literal["min", "max"] = "min" + ) -> Result: + """ + Re-solve an existing native solver model and return a Result. + + This enables iterative solve workflows: solve once via + ``Model.solve()``, modify the native ``solver_model`` (e.g. change + objective coefficients), then call ``resolve()`` followed by + ``Model.apply_result()``. + + Parameters + ---------- + solver_model : Any + The native solver object (e.g. ``highspy.Highs``, ``gurobipy.Model``). + sense : {{"min", "max"}} + Optimization sense. + + Returns + ------- + Result + """ + msg = f"{type(self).__name__} does not support resolve()" + raise NotImplementedError(msg) + @property def solver_name(self) -> SolverName: return SolverName[self.__class__.__name__] @@ -920,6 +945,11 @@ def solve_problem_from_file( sense=read_sense_from_problem_file(problem_fn), ) + def resolve( + self, solver_model: highspy.Highs, sense: Literal["min", "max"] = "min" + ) -> Result: + return self._solve(solver_model, io_api="direct", sense=sense) + def _set_solver_params( self, highs_solver: highspy.Highs, @@ -1162,6 +1192,19 @@ def solve_problem_from_file( sense=sense, ) + def resolve( + self, solver_model: gurobipy.Model, sense: Literal["min", "max"] = "min" + ) -> Result: + return self._solve( + solver_model, + solution_fn=None, + log_fn=None, + warmstart_fn=None, + basis_fn=None, + io_api="direct", + sense=sense, + ) + def _solve( self, m: gurobipy.Model, @@ -2094,6 +2137,19 @@ def solve_problem_from_file( sense=sense, ) + def resolve( + self, solver_model: mosek.Task, sense: Literal["min", "max"] = "min" + ) -> Result: + return self._solve( + solver_model, + solution_fn=None, + log_fn=None, + warmstart_fn=None, + basis_fn=None, + io_api="direct", + sense=sense, + ) + def _solve( self, m: mosek.Task, diff --git a/test/test_resolve.py b/test/test_resolve.py new file mode 100644 index 00000000..913b5249 --- /dev/null +++ b/test/test_resolve.py @@ -0,0 +1,111 @@ +from __future__ import annotations + +from pathlib import Path +from typing import Any, Literal, cast + +import numpy as np +import pytest + +from linopy import Model +from linopy.constants import Result, Solution, Status, TerminationCondition +from linopy.solvers import Highs, SolverName + + +@pytest.fixture +def model() -> Model: + m = Model() + x = m.add_variables(lower=0, upper=10, name="x") + y = m.add_variables(lower=0, upper=10, name="y") + m.add_constraints(x + y >= 8, name="c1") + m.add_constraints(x + y <= 15, name="c2") + m.objective = 2 * x + 3 * y + return m + + +def test_apply_result_basic(model: Model) -> None: + model.solve(solver_name="highs") + + status = Status.from_termination_condition(TerminationCondition.optimal) + solution = Solution(objective=42.0) + result = Result(status=status, solution=solution) + + s, tc = model.apply_result(result, solver_name="test") + assert s == "ok" + assert model.objective.value == 42.0 + assert model.solver_name == "test" + + +def test_apply_result_infeasible(model: Model) -> None: + status = Status.from_termination_condition(TerminationCondition.infeasible) + result = Result(status=status, solution=Solution(objective=np.nan)) + + s, tc = model.apply_result(result) + assert s == "warning" + assert tc == "infeasible" + assert model.solver_name == "unknown" + + +def test_apply_result_none_solution(model: Model) -> None: + status = Status.from_termination_condition(TerminationCondition.infeasible) + result = Result(status=status, solution=None) + + s, tc = model.apply_result(result) + assert s == "warning" + assert tc == "infeasible" + + +def test_resolve_highs(model: Model) -> None: + pytest.importorskip("highspy") + model.solve(solver_name="highs") + obj1 = model.objective.value + + h = model.solver_model + n_cols = h.getNumCol() + new_costs = np.array([10.0, 1.0], dtype=float) + h.changeColsCost(n_cols, np.arange(n_cols, dtype=np.int32), new_costs) + + solver = Highs() + sense = cast(Literal["min", "max"], model.sense) + result = solver.resolve(h, sense=sense) + model.apply_result(result, solver_name="highs") + + assert model.status == "ok" + assert model.objective.value != obj1 + assert model.solution["y"].values.item() >= model.solution["x"].values.item() + + +def test_resolve_not_supported() -> None: + from linopy.solvers import Solver + + class DummySolver(Solver[None]): + def solve_problem_from_model( + self, + model: Model, + solution_fn: Path | None = None, + log_fn: Path | None = None, + warmstart_fn: Path | None = None, + basis_fn: Path | None = None, + env: Any = None, + explicit_coordinate_names: bool = False, + ) -> Result: + raise NotImplementedError + + def solve_problem_from_file( + self, + problem_fn: Path, + solution_fn: Path | None = None, + log_fn: Path | None = None, + warmstart_fn: Path | None = None, + basis_fn: Path | None = None, + env: Any = None, + ) -> Result: + raise NotImplementedError + + @property + def solver_name(self) -> SolverName: + return SolverName.Highs + + solver = DummySolver.__new__(DummySolver) + solver.solver_options = {} + with pytest.raises(NotImplementedError, match="does not support resolve"): + solver.resolve(None)