diff --git a/CLAUDE.md b/CLAUDE.md
index 1f696a0b..e5e3c182 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -7,22 +7,22 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 ### Running Tests
 ```bash
 # Run all tests (excluding GPU tests by default)
-pytest
+uv run --extra dev --extra solvers pytest
 
 # Run tests with coverage
-pytest --cov=./ --cov-report=xml linopy --doctest-modules test
+uv run --extra dev --extra solvers pytest --cov=./ --cov-report=xml linopy --doctest-modules test
 
 # Run a specific test file
-pytest test/test_model.py
+uv run --extra dev --extra solvers pytest test/test_model.py
 
 # Run a specific test function
-pytest test/test_model.py::test_model_creation
+uv run --extra dev --extra solvers pytest test/test_model.py::test_model_creation
 
 # Run GPU tests (requires GPU hardware and cuPDLPx installation)
-pytest --run-gpu
+uv run --extra dev --extra solvers pytest --run-gpu
 
 # Run only GPU tests
-pytest -m gpu --run-gpu
+uv run --extra dev --extra solvers pytest -m gpu --run-gpu
 ```
 
 **GPU Testing**: Tests that require GPU hardware (e.g., cuPDLPx solver) are automatically skipped by default since CI machines typically don't have GPUs. To run GPU tests locally, use the `--run-gpu` flag. The tests are automatically marked with `@pytest.mark.gpu` based on solver capabilities.
@@ -30,17 +30,17 @@ pytest -m gpu --run-gpu
 ### Linting and Type Checking
 ```bash
 # Run linter (ruff)
-ruff check .
-ruff check --fix .  # Auto-fix issues
+uv run --extra dev ruff check .
+uv run --extra dev ruff check --fix .  # Auto-fix issues
 
 # Run formatter
-ruff format .
+uv run --extra dev ruff format .
 
 # Run type checker
-mypy .
+uv run --extra dev mypy .
 
 # Run all pre-commit hooks
-pre-commit run --all-files
+uv run --extra dev pre-commit run --all-files
 ```
 
 ### Development Setup
diff --git a/dev-scripts/piecewise-feasibility-tests-walkthrough.ipynb b/dev-scripts/piecewise-feasibility-tests-walkthrough.ipynb
deleted file mode 100644
index b2fdaf7c..00000000
--- a/dev-scripts/piecewise-feasibility-tests-walkthrough.ipynb
+++ /dev/null
@@ -1,348 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# `test_piecewise_feasibility.py` — visual walkthrough\n",
-    "\n",
-    "**Purpose:** document what each test class in `test/test_piecewise_feasibility.py` actually probes, with pictures. Intended as review aid for the PR — **not** merged into master.\n",
-    "\n",
-    "The test file stress-tests the claim that `add_piecewise_formulation(sign=\"<=\"/\">=\")` yields the **same feasible region** for `(x, y)` regardless of which method (`lp` / `sos2` / `incremental`) dispatches the formulation, on curves where all three are applicable.\n",
-    "\n",
-    "Four test classes:\n",
-    "\n",
-    "| class | what it probes | scope |\n",
-    "|---|---|---|\n",
-    "| `TestRotatedObjective` | support-function equivalence — 16 rotation directions | the strong test |\n",
-    "| `TestDomainBoundary` | `x` outside `[x_min, x_max]` is infeasible | LP explicit vs SOS2 implicit |\n",
-    "| `TestPointwiseInfeasibility` | `y` just past `f(x)` is infeasible | targeted sanity check |\n",
-    "| `TestNVariableInequality` | 3-variable: first tuple bounded, rest equality | SOS2 vs incremental only |\n",
-    "\n",
-    "Below: one visualization per class.\n",
-    "\n",
-    "*Run this notebook from the repository root so that `from test.test_piecewise_feasibility import ...` resolves.*"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-04-23T08:00:47.376525Z",
-     "start_time": "2026-04-23T08:00:46.142492Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
-    "\n",
-    "from test.test_piecewise_feasibility import (\n",
-    "    CURVES,\n",
-    "    Y_HI,\n",
-    "    Y_LO,\n",
-    "    Curve,\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Shared primitive: draw the curve and its feasible region"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-04-23T08:00:47.384959Z",
-     "start_time": "2026-04-23T08:00:47.381361Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "def draw_curve_and_region(ax, curve: Curve, *, shade: bool = True) -> None:\n",
-    "    \"\"\"Plot breakpoints + shade the feasible region (hypograph or epigraph).\"\"\"\n",
-    "    xs = np.array(curve.x_pts)\n",
-    "    ys = np.array(curve.y_pts)\n",
-    "    ax.plot(xs, ys, \"o-\", color=\"C0\", lw=2, label=\"breakpoints\")\n",
-    "\n",
-    "    if shade:\n",
-    "        if curve.sign == \"<=\":\n",
-    "            ax.fill_between(\n",
-    "                xs,\n",
-    "                np.full_like(ys, Y_LO),\n",
-    "                ys,\n",
-    "                alpha=0.15,\n",
-    "                color=\"C0\",\n",
-    "                label=f\"feasible: y {curve.sign} f(x)\",\n",
-    "            )\n",
-    "        else:\n",
-    "            ax.fill_between(\n",
-    "                xs,\n",
-    "                ys,\n",
-    "                np.full_like(ys, Y_HI),\n",
-    "                alpha=0.15,\n",
-    "                color=\"C0\",\n",
-    "                label=f\"feasible: y {curve.sign} f(x)\",\n",
-    "            )\n",
-    "\n",
-    "    pad_x = 0.15 * (xs.max() - xs.min())\n",
-    "    pad_y = 0.15 * (ys.max() - ys.min()) + 1\n",
-    "    ax.set_xlim(xs.min() - pad_x, xs.max() + pad_x)\n",
-    "    ax.set_ylim(ys.min() - pad_y, ys.max() + pad_y)\n",
-    "    ax.set_xlabel(\"x\")\n",
-    "    ax.set_ylabel(\"y\")\n",
-    "    ax.grid(alpha=0.3)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## `TestRotatedObjective` — the strong test\n",
-    "\n",
-    "For every direction `(α, β)` on the unit circle, minimize `α·x + β·y` under the PWL. The answer is the **support function** of the feasible region in direction `(α, β)` — and for a convex region, the support function uniquely determines the region. If LP and SOS2/incremental give the same support-function value for 16 directions, their feasible regions are identical.\n",
-    "\n",
-    "Each red dot below is the extreme point the solver lands at for one direction. The arrows show the objective-push direction. A failure would manifest as one method's dot landing at a different vertex than the oracle's."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-04-23T08:00:47.829483Z",
-     "start_time": "2026-04-23T08:00:47.388542Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "def panel_rotated_objective(ax, curve: Curve, n_dirs: int = 16) -> None:\n",
-    "    draw_curve_and_region(ax, curve)\n",
-    "    xs, ys = np.array(curve.x_pts), np.array(curve.y_pts)\n",
-    "    cx = 0.5 * (xs.min() + xs.max())\n",
-    "    cy = 0.5 * (ys.min() + ys.max())\n",
-    "    arrow_len = 0.25 * min(xs.max() - xs.min(), (ys.max() - ys.min()) + 5)\n",
-    "\n",
-    "    for i in range(n_dirs):\n",
-    "        theta = 2 * np.pi * i / n_dirs\n",
-    "        alpha, beta = np.cos(theta), np.sin(theta)\n",
-    "        ax.annotate(\n",
-    "            \"\",\n",
-    "            xytext=(cx, cy),\n",
-    "            xy=(cx + arrow_len * alpha, cy + arrow_len * beta),\n",
-    "            arrowprops=dict(arrowstyle=\"->\", color=\"C3\", alpha=0.4, lw=1),\n",
-    "        )\n",
-    "        # Oracle extreme point in this direction\n",
-    "        verts = curve.vertices()\n",
-    "        extreme = min(verts, key=lambda v: alpha * v[0] + beta * v[1])\n",
-    "        ax.plot(*extreme, \"o\", color=\"C3\", ms=4, alpha=0.7)\n",
-    "\n",
-    "    ax.plot([], [], \"o\", color=\"C3\", alpha=0.7, label=f\"{n_dirs} extreme points\")\n",
-    "    ax.legend(loc=\"upper left\", fontsize=8)\n",
-    "    ax.set_title(f\"{curve.name}  (sign={curve.sign})\")\n",
-    "\n",
-    "\n",
-    "fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))\n",
-    "panel_rotated_objective(axes[0], CURVES[0])  # concave-smooth\n",
-    "panel_rotated_objective(axes[1], CURVES[2])  # convex-steep\n",
-    "panel_rotated_objective(axes[2], CURVES[5])  # two-segment\n",
-    "fig.suptitle(\n",
-    "    \"TestRotatedObjective — support function sampled at 16 directions\", fontsize=12\n",
-    ")\n",
-    "plt.tight_layout();"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Notice the **dots cluster at the curve breakpoints** (top edges) and at the **bottom corners** `(x_min, Y_LO)`, `(x_max, Y_LO)`. That's because the feasible region is a polygon: linear objectives always attain their optimum at a vertex.\n",
-    "\n",
-    "The 288 pytest items (6 curves × 3 methods × 16 directions) check that all three methods land at the same extreme point for every direction."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## `TestDomainBoundary` — enforce `x ∈ [x_min, x_max]`\n",
-    "\n",
-    "LP enforces this with an explicit constraint; SOS2/incremental enforce it implicitly via `sum(λ) = 1`. Two different implementations of the same bound — worth a direct probe."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-04-23T08:00:48.103641Z",
-     "start_time": "2026-04-23T08:00:47.835275Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "def panel_domain_boundary(ax, curve: Curve) -> None:\n",
-    "    draw_curve_and_region(ax, curve)\n",
-    "    xs = np.array(curve.x_pts)\n",
-    "    y_span = ax.get_ylim()\n",
-    "    ax.axvline(xs[0], color=\"C2\", lw=1.5, label=f\"x_min={xs[0]}\")\n",
-    "    ax.axvline(xs[-1], color=\"C2\", lw=1.5, label=f\"x_max={xs[-1]}\")\n",
-    "    ax.axvline(xs[0] - 1, color=\"C3\", lw=1.5, ls=\"--\")\n",
-    "    ax.axvline(xs[-1] + 1, color=\"C3\", lw=1.5, ls=\"--\")\n",
-    "    yy = y_span[1] - 0.12 * (y_span[1] - y_span[0])\n",
-    "    ax.text(\n",
-    "        xs[0] - 1, yy, \"INFEASIBLE\\n(x < x_min)\", ha=\"center\", fontsize=8, color=\"C3\"\n",
-    "    )\n",
-    "    ax.text(\n",
-    "        xs[-1] + 1, yy, \"INFEASIBLE\\n(x > x_max)\", ha=\"center\", fontsize=8, color=\"C3\"\n",
-    "    )\n",
-    "    ax.legend(loc=\"lower center\", fontsize=7)\n",
-    "    ax.set_title(f\"{curve.name} — domain probe\")\n",
-    "\n",
-    "\n",
-    "fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))\n",
-    "panel_domain_boundary(axes[0], CURVES[0])  # concave-smooth\n",
-    "panel_domain_boundary(axes[1], CURVES[1])  # concave-shifted (negative domain)\n",
-    "panel_domain_boundary(axes[2], CURVES[5])  # two-segment\n",
-    "fig.suptitle(\n",
-    "    \"TestDomainBoundary — x outside the breakpoint range is infeasible\", fontsize=12\n",
-    ")\n",
-    "plt.tight_layout();"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## `TestPointwiseInfeasibility` — y just past the curve\n",
-    "\n",
-    "Rotated objectives probe *extremes*; this test specifically nudges `y` past `f(x)` by a small margin (`0.01`) and asserts infeasibility. Catches NaN-mask or off-by-one-segment bugs that might accidentally allow slack."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-04-23T08:00:48.366674Z",
-     "start_time": "2026-04-23T08:00:48.112127Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "def panel_pointwise(ax, curve: Curve) -> None:\n",
-    "    draw_curve_and_region(ax, curve)\n",
-    "    xs = np.array(curve.x_pts)\n",
-    "    x_mid = 0.5 * (xs[0] + xs[-1])\n",
-    "    fx = curve.f(x_mid)\n",
-    "    y_bad = fx + 0.01 if curve.sign == \"<=\" else fx - 0.01\n",
-    "    ax.plot(x_mid, fx, \"o\", color=\"C2\", ms=9, label=f\"on curve: f({x_mid:g})={fx:g}\")\n",
-    "    ax.plot(\n",
-    "        x_mid, y_bad, \"x\", color=\"C3\", ms=14, mew=3, label=f\"infeasible: y={y_bad:g}\"\n",
-    "    )\n",
-    "    ax.legend(loc=\"lower right\", fontsize=7)\n",
-    "    ax.set_title(f\"{curve.name} — nudge past f(x)\")\n",
-    "\n",
-    "\n",
-    "fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))\n",
-    "panel_pointwise(axes[0], CURVES[0])  # concave-smooth, sign=\"<=\"\n",
-    "panel_pointwise(axes[1], CURVES[2])  # convex-steep, sign=\">=\"\n",
-    "panel_pointwise(axes[2], CURVES[4])  # linear-gte\n",
-    "fig.suptitle(\n",
-    "    \"TestPointwiseInfeasibility — y past the curve by 0.01 in the sign direction\",\n",
-    "    fontsize=12,\n",
-    ")\n",
-    "plt.tight_layout();"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## `TestNVariableInequality` — 3-variable sign split\n",
-    "\n",
-    "With three tuples `(fuel, power, heat)` and `sign=\"<=\"`:\n",
-    "- `fuel` (the **first** tuple) is **bounded above** by its interpolated value,\n",
-    "- `power` and `heat` (remaining tuples) are **forced to equality** — pinned on the curve.\n",
-    "\n",
-    "LP doesn't support N > 2 tuples, so this class compares SOS2 vs incremental only. The 3D plot shows the CHP curve and the 7 test points (one per `power_fix`) that both methods must agree on."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-04-23T08:00:48.489668Z",
-     "start_time": "2026-04-23T08:00:48.371526Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "bp = {\n",
-    "    \"power\": np.array([0, 30, 60, 100]),\n",
-    "    \"fuel\": np.array([0, 40, 85, 160]),\n",
-    "    \"heat\": np.array([0, 25, 55, 95]),\n",
-    "}\n",
-    "\n",
-    "fig = plt.figure(figsize=(9, 6.5))\n",
-    "ax = fig.add_subplot(projection=\"3d\")\n",
-    "ax.plot(\n",
-    "    bp[\"power\"], bp[\"fuel\"], bp[\"heat\"], \"o-\", color=\"C0\", lw=2, label=\"CHP breakpoints\"\n",
-    ")\n",
-    "\n",
-    "for p in [0, 15, 30, 45, 60, 80, 100]:\n",
-    "    f = np.interp(p, bp[\"power\"], bp[\"fuel\"])\n",
-    "    h = np.interp(p, bp[\"power\"], bp[\"heat\"])\n",
-    "    ax.plot([p], [f], [h], \"o\", color=\"C3\", ms=7)\n",
-    "    # drop to base plane\n",
-    "    ax.plot([p, p], [f, 0], [h, h], color=\"C3\", alpha=0.3, lw=0.8)\n",
-    "\n",
-    "ax.set_xlabel(\"power\")\n",
-    "ax.set_ylabel(\"fuel\")\n",
-    "ax.set_zlabel(\"heat\")\n",
-    "ax.plot(\n",
-    "    [],\n",
-    "    [],\n",
-    "    \"o\",\n",
-    "    color=\"C3\",\n",
-    "    label=\"7 test points — power pinned,\\nfuel at upper bound, heat on curve\",\n",
-    ")\n",
-    "ax.set_title('TestNVariableInequality — CHP curve (sign=\"<=\")')\n",
-    "ax.legend(loc=\"upper left\", fontsize=8);"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## What a failing test would tell you\n",
-    "\n",
-    "- **Rotated objective fails**: the methods disagree on the feasible region in some direction. The failure message includes the attained `(x, y)` point — you'd see which extreme point one method landed at that the others didn't.\n",
-    "- **Domain boundary fails**: one method lets `x` escape `[x_min, x_max]`. LP path most likely: the domain-bound constraint was dropped. SOS2 path: the `sum(λ) = 1` constraint was weakened.\n",
-    "- **Pointwise infeasibility fails**: one method accepts a point past the curve. Most often a NaN-mask bug in per-entity formulations, or a wrong segment getting picked.\n",
-    "- **N-variable fails**: the sign split went wrong — either an input leaked into the signed link or the first-tuple convention is misrouting.\n",
-    "\n",
-    "All 356 pytest items are currently green at `TOL = 1e-5`."
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "name": "python",
-   "version": "3.13.2"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 75c316a0..d399afe3 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -4,6 +4,12 @@ Release Notes
 Upcoming Version
 ----------------
 
+* Solver refactor: solver state now lives on a stateful ``Solver`` instance attached to ``Model.solver``. ``Model.solver_model`` and ``Model.solver_name`` become read-only properties delegating to ``model.solver`` (assigning anything other than ``None`` raises; setting ``None`` closes the solver). ``Model.solver_name`` may be ``None`` before a solve. The latter two properties may be deprecated in future versions.
+* Advanced solve workflow: construct via ``Solver.from_name(name, model, io_api=..., options=...)`` (or ``SolverClass.from_model(model, ...)``), then call ``solver.solve()`` to run and obtain a ``Result``, and ``model.apply_result(result)`` to write the solution back to the model. ``Solver`` is now a dataclass; subclasses no longer need ``__init__`` overrides. The previous two-step ``Model.prepare_solver`` / ``Model.run_solver`` API has been removed (it was added in the same upcoming release and not yet shipped).
+* Solver capabilities are declared as ``features: frozenset[SolverFeature]`` ClassVars on each ``Solver`` subclass; use ``Solver.supports(feature)``. ``SolverFeature`` is now exported from ``linopy`` (and from ``linopy.solvers``); ``linopy.solver_capabilities`` remains as a back-compat shim with a lazy ``SOLVER_REGISTRY`` mapping.
+* ``Result`` gains ``solver_name`` and ``report: SolverReport | None`` (runtime, MIP gap, dual bound, iteration counts) and prints them in ``__repr__``. CBC, HiGHS, Gurobi, Knitro, and cuPDLPx populate ``report``; when possible also populate the MIP ``dual_bound``.
+* ``Solution.primal`` and ``Solution.dual`` are now ``np.ndarray`` lookup arrays keyed by integer model labels (length = ``max_label + 1``, ``NaN`` for unfilled positions); previously ``pd.Series`` keyed by variable/constraint name. Each solver is responsible for emitting the label-indexed form — direct-API solvers via cached ``_vlabels``/``_clabels`` populated at ``_build_direct`` time, file-based solvers via the shared ``_solution_from_names`` helper which parses linopy labels from solver-side names. Fixes solution mapping for file-based solvers that may iterate variables in a different order than linopy's build order or drop unused variables entirely. Internal — code that introspected these must use ``np.ndarray`` semantics. Solution mapping reads labels from a cached ``ConstraintLabelIndex`` on ``Model.constraints`` and no longer triggers a constraint-matrix rebuild.
+* Per-solver translation helpers are unified under ``Solver._build_solver_model``. The module-level ``to_gurobipy`` / ``to_highspy`` / ``to_mosek`` / ``to_cupdlpx`` (and their ``Model`` bindings) are kept as thin wrappers that route through ``Solver.from_model(model, io_api="direct")`` and return the native solver model.
 * Increase speed of direct solver communication (~10x) in conversion functions like `to_highspy` through faster matrix creation (see below), leading to significant overall speed-up when setting `io_api="direct"`.
 * Add ``CSRConstraint``, a memory-efficient immutable constraint representation backed by scipy CSR sparse matrices. Provides up to 90% memory savings for constraints with many terms and 30–120x faster matrix generation for direct solver APIs.
   - Add ``freeze_constraints`` parameter to ``Model`` for globally storing constraints in CSR format on ``add_constraints``.
diff --git a/linopy/__init__.py b/linopy/__init__.py
index df07cc81..cf2e7832 100644
--- a/linopy/__init__.py
+++ b/linopy/__init__.py
@@ -39,6 +39,7 @@
     tangent_lines,
 )
 from linopy.remote import RemoteHandler
+from linopy.solvers import SolverFeature
 
 try:
     from linopy.remote import OetcCredentials, OetcHandler, OetcSettings  # noqa: F401
@@ -63,6 +64,7 @@
     "QuadraticExpression",
     "RemoteHandler",
     "Slopes",
+    "SolverFeature",
     "Variable",
     "Variables",
     "align",
diff --git a/linopy/common.py b/linopy/common.py
index 162fcdfe..e9a38d29 100644
--- a/linopy/common.py
+++ b/linopy/common.py
@@ -999,6 +999,43 @@ def invalidate(self) -> None:
         self.__dict__.pop("label_to_pos", None)
 
 
+class ConstraintLabelIndex:
+    """
+    Index for O(1) mapping between constraint labels and dense positions.
+
+    Mirrors VariableLabelIndex on the constraint side, but without building
+    the full constraint matrix — only labels and the row mask are computed.
+    """
+
+    def __init__(self, constraints: Any) -> None:
+        self._constraints = constraints
+
+    @cached_property
+    def clabels(self) -> np.ndarray:
+        """Active constraint labels in build order, shape (n_active_cons,)."""
+        label_lists = [c.active_labels() for c in self._constraints.data.values()]
+        return (
+            np.concatenate(label_lists) if label_lists else np.array([], dtype=np.intp)
+        )
+
+    @cached_property
+    def label_to_pos(self) -> np.ndarray:
+        """Mapping from constraint label to dense position, shape (_cCounter,)."""
+        clabels = self.clabels
+        n = self._constraints.model._cCounter
+        label_to_pos = np.full(n, -1, dtype=np.intp)
+        label_to_pos[clabels] = np.arange(len(clabels), dtype=np.intp)
+        return label_to_pos
+
+    @property
+    def n_active_cons(self) -> int:
+        return len(self.clabels)
+
+    def invalidate(self) -> None:
+        self.__dict__.pop("clabels", None)
+        self.__dict__.pop("label_to_pos", None)
+
+
 def get_label_position(
     obj: Any,
     values: int | np.ndarray,
@@ -1509,49 +1546,34 @@ def is_constant(x: SideLike) -> bool:
     )
 
 
-def series_to_lookup_array(s: pd.Series) -> np.ndarray:
+def values_to_lookup_array(
+    values: np.ndarray, labels: np.ndarray, size: int | None = None
+) -> np.ndarray:
     """
-    Convert an integer-indexed Series to a dense numpy lookup array.
+    Build a dense NaN-padded lookup array from values and integer labels.
 
-    Non-negative indices are placed at their corresponding positions;
-    negative indices are ignored. Gaps are filled with NaN.
+    Non-negative labels are placed at their corresponding positions; negative
+    labels are skipped. Gaps are filled with NaN.
 
     Parameters
     ----------
-    s : pd.Series
-        Series with an integer index.
+    values : np.ndarray
+        Values to place into the lookup array.
+    labels : np.ndarray
+        Integer labels giving the target position for each value.
+    size : int, optional
+        Length of the returned array. Defaults to ``max(labels) + 1`` if any
+        non-negative label is present, otherwise 0.
 
     Returns
     -------
     np.ndarray
-        Dense array of length ``max(index) + 1``.
-    """
-    max_idx = max(int(s.index.max()), 0)
-    arr = np.full(max_idx + 1, nan)
-    mask = s.index >= 0
-    arr[s.index[mask]] = s.values[mask]
+        Dense float lookup array.
+    """
+    labels = np.asarray(labels, dtype=int)
+    mask = labels >= 0
+    if size is None:
+        size = int(labels[mask].max()) + 1 if mask.any() else 0
+    arr = np.full(size, nan, dtype=float)
+    arr[labels[mask]] = values[mask]
     return arr
-
-
-def lookup_vals(arr: np.ndarray, idx: np.ndarray) -> np.ndarray:
-    """
-    Look up values from a dense array by integer labels.
-
-    Negative labels and labels beyond the array length map to NaN.
-
-    Parameters
-    ----------
-    arr : np.ndarray
-        Dense lookup array (e.g. from :func:`series_to_lookup_array`).
-    idx : np.ndarray
-        Integer label indices.
-
-    Returns
-    -------
-    np.ndarray
-        Array of looked-up values with the same shape as *idx*.
-    """
-    valid = (idx >= 0) & (idx < len(arr))
-    vals = np.full(idx.shape, nan)
-    vals[valid] = arr[idx[valid]]
-    return vals
diff --git a/linopy/constants.py b/linopy/constants.py
index 5cc98ce2..86af18ce 100644
--- a/linopy/constants.py
+++ b/linopy/constants.py
@@ -9,7 +9,6 @@
 from typing import Any, Literal, TypeAlias, Union, get_args
 
 import numpy as np
-import pandas as pd
 
 logger = logging.getLogger(__name__)
 
@@ -261,21 +260,35 @@ def is_ok(self) -> bool:
         return self.status == SolverStatus.ok
 
 
-def _pd_series_float() -> pd.Series:
-    return pd.Series(dtype=float)
-
-
 @dataclass
 class Solution:
     """
     Solution returned by the solver.
+
+    ``primal`` and ``dual`` are dense float arrays indexed by linopy label:
+    ``primal[label]`` is the value for variable ``label``, with ``NaN`` where
+    no value is available (masked labels, vars dropped by the solver, etc.).
+    Each solver is responsible for emitting arrays in this label-indexed form.
     """
 
-    primal: pd.Series = field(default_factory=_pd_series_float)
-    dual: pd.Series = field(default_factory=_pd_series_float)
+    primal: np.ndarray = field(default_factory=lambda: np.array([], dtype=float))
+    dual: np.ndarray = field(default_factory=lambda: np.array([], dtype=float))
     objective: float = field(default=np.nan)
 
 
+@dataclass
+class SolverReport:
+    """
+    Solver-reported performance metrics.
+    """
+
+    runtime: float | None = None
+    mip_gap: float | None = None
+    dual_bound: float | None = None
+    barrier_iterations: int | None = None
+    simplex_iterations: int | None = None
+
+
 @dataclass
 class Result:
     """
@@ -285,6 +298,8 @@ class Result:
     status: Status
     solution: Solution | None = None
     solver_model: Any = None
+    solver_name: str = ""
+    report: SolverReport | None = None
 
     def __repr__(self) -> str:
         solver_model_string = (
@@ -297,10 +312,21 @@ def __repr__(self) -> str:
             )
         else:
             solution_string = "Solution: None\n"
+        solver_name_string = f"Solver: {self.solver_name}\n" if self.solver_name else ""
+        report_string = ""
+        if self.report is not None:
+            if self.report.runtime is not None:
+                report_string += f"Runtime: {self.report.runtime:.2f}s\n"
+            if self.report.mip_gap is not None:
+                report_string += f"MIP gap: {self.report.mip_gap:.2e}\n"
+            if self.report.dual_bound is not None:
+                report_string += f"Dual bound: {self.report.dual_bound:.2e}\n"
         return (
             f"Status: {self.status.status.value}\n"
             f"Termination condition: {self.status.termination_condition.value}\n"
             + solution_string
+            + solver_name_string
+            + report_string
             + f"Solver model: {solver_model_string}\n"
             f"Solver message: {self.status.legacy_status}"
         )
diff --git a/linopy/constraints.py b/linopy/constraints.py
index 6aab4902..b74dee5c 100644
--- a/linopy/constraints.py
+++ b/linopy/constraints.py
@@ -32,6 +32,7 @@
 
 from linopy import expressions, variables
 from linopy.common import (
+    ConstraintLabelIndex,
     LabelPositionIndex,
     LocIndexer,
     VariableLabelIndex,
@@ -143,6 +144,11 @@ def is_assigned(self) -> bool:
     def labels(self) -> DataArray:
         """Get the labels DataArray."""
 
+    @property
+    @abstractmethod
+    def range(self) -> tuple[int, int]:
+        """Return the label range of the constraint."""
+
     @property
     @abstractmethod
     def coeffs(self) -> DataArray:
@@ -212,6 +218,10 @@ def to_matrix_with_rhs(
         the RHS/sense vectors are needed.
         """
 
+    @abstractmethod
+    def active_labels(self) -> np.ndarray:
+        """Active constraint labels in build order, without building the CSR."""
+
     def __getitem__(
         self, selector: str | int | slice | list | tuple | dict
     ) -> Constraint:
@@ -561,6 +571,10 @@ def attrs(self) -> dict[str, Any]:
             d["label_range"] = (self._cindex, self._cindex + self.full_size)
         return d
 
+    @property
+    def coords(self) -> DatasetCoordinates:
+        return Dataset(coords={c.name: c for c in self._coords}).coords
+
     @property
     def dims(self) -> Frozen[Hashable, int]:
         d: dict[Hashable, int] = {c.name: len(c) for c in self._coords}
@@ -865,6 +879,9 @@ def to_matrix_with_rhs(
             sense = np.array([s[0] for s in self._sign])
         return self._csr, self._con_labels, self._rhs, sense
 
+    def active_labels(self) -> np.ndarray:
+        return self._con_labels
+
     def sanitize_zeros(self) -> CSRConstraint:
         """Remove terms with zero or near-zero coefficients (mutates in-place)."""
         self._csr.data[np.abs(self._csr.data) <= 1e-10] = 0
@@ -1222,6 +1239,18 @@ def to_matrix(
         csr.sum_duplicates()
         return csr, con_labels
 
+    def active_labels(self) -> np.ndarray:
+        labels_flat = self.labels.values.ravel()
+        vars_vals = self.vars.values
+        n_rows = len(labels_flat)
+        vars_2d = (
+            vars_vals.reshape(n_rows, -1)
+            if n_rows > 0
+            else vars_vals.reshape(0, max(1, vars_vals.size))
+        )
+        row_mask = (labels_flat != -1) & (vars_2d != -1).any(axis=1)
+        return labels_flat[row_mask]
+
     def to_matrix_with_rhs(
         self, label_index: VariableLabelIndex
     ) -> tuple[scipy.sparse.csr_array, np.ndarray, np.ndarray, np.ndarray]:
@@ -1427,6 +1456,7 @@ class Constraints:
     data: dict[str, ConstraintBase]
     model: Model
     _label_position_index: LabelPositionIndex | None = None
+    _constraint_label_index: ConstraintLabelIndex | None = None
 
     dataset_attrs = ["labels", "coeffs", "vars", "sign", "rhs"]
     dataset_names = [
@@ -1548,6 +1578,15 @@ def _invalidate_label_position_index(self) -> None:
         """Invalidate the label position index cache."""
         if self._label_position_index is not None:
             self._label_position_index.invalidate()
+        if self._constraint_label_index is not None:
+            self._constraint_label_index.invalidate()
+
+    @property
+    def label_index(self) -> ConstraintLabelIndex:
+        """Index for O(1) label->position mapping and compact clabels array."""
+        if self._constraint_label_index is None:
+            self._constraint_label_index = ConstraintLabelIndex(self)
+        return self._constraint_label_index
 
     @property
     def labels(self) -> Dataset:
diff --git a/linopy/io.py b/linopy/io.py
index 6dc1c9c9..36d7abb3 100644
--- a/linopy/io.py
+++ b/linopy/io.py
@@ -9,7 +9,6 @@
 import logging
 import shutil
 import time
-import warnings
 from collections.abc import Callable, Iterable
 from io import BufferedWriter
 from pathlib import Path
@@ -20,8 +19,6 @@
 import pandas as pd
 import polars as pl
 import xarray as xr
-from numpy import ones_like, zeros_like
-from scipy.sparse import tril, triu
 from tqdm import tqdm
 
 from linopy import solvers
@@ -625,7 +622,9 @@ def to_file(
 
         # Use very fast highspy implementation
         # Might be replaced by custom writer, however needs C/Rust bindings for performance
-        h = m.to_highspy(explicit_coordinate_names=explicit_coordinate_names)
+        h = solvers.Highs._build_solver_model(
+            m, explicit_coordinate_names=explicit_coordinate_names
+        )
         h.writeModel(str(fn))
     else:
         raise ValueError(
@@ -641,125 +640,17 @@ def to_mosek(
     explicit_coordinate_names: bool = False,
     set_names: bool = True,
 ) -> Any:
-    """
-    Export model to MOSEK.
-
-    Export the model directly to MOSEK without writing files.
-
-    Parameters
-    ----------
-    m : linopy.Model
-    task : empty MOSEK task
-    explicit_coordinate_names : bool, optional
-        Whether to use explicit coordinate names. Default is False.
-    set_names : bool, optional
-        Whether to set variable and constraint names. Default is True.
-        Setting to False can significantly speed up model export.
-
-    Returns
-    -------
-    task : MOSEK Task object
-    """
-    if m.variables.sos:
-        raise NotImplementedError("SOS constraints are not supported by MOSEK.")
-
-    if m.variables.semi_continuous:
-        raise NotImplementedError(
-            "Semi-continuous variables are not supported by MOSEK. "
-            "Use a solver that supports them (gurobi, cplex, highs)."
-        )
-
+    """Build the MOSEK task for `m`."""
     import mosek
 
     if task is None:
         task = mosek.Task()
-
-    task.appendvars(m.nvars)
-    task.appendcons(m.ncons)
-
-    M = m.matrices
-
-    if set_names:
-        print_variables, print_constraints = get_printers_scalar(
-            m, explicit_coordinate_names=explicit_coordinate_names
-        )
-        labels = print_variables(M.vlabels)
-        task.generatevarnames(
-            np.arange(0, len(labels)), "%0", [len(labels)], None, [0], labels
-        )
-
-    ## Variables
-
-    # MOSEK uses bound keys (free, bounded below or above, ranged and fixed)
-    # plus bound values (lower and upper), and it is considered an error to
-    # input an infinite value for a finite bound.
-    # bkx and bkc define the boundkeys based on upper and lower bound, and blx,
-    # bux, blc and buc define the finite bounds. The numerical value of a bound
-    # indicated to be infinite by the bound key is ignored by MOSEK.
-    bkx = [
-        (
-            (
-                (mosek.boundkey.ra if lb < ub else mosek.boundkey.fx)
-                if ub < np.inf
-                else mosek.boundkey.lo
-            )
-            if (lb > -np.inf)
-            else (mosek.boundkey.up if (ub < np.inf) else mosek.boundkey.fr)
-        )
-        for (lb, ub) in zip(M.lb, M.ub)
-    ]
-    blx = [b if b > -np.inf else 0.0 for b in M.lb]
-    bux = [b if b < np.inf else 0.0 for b in M.ub]
-    task.putvarboundslice(0, m.nvars, bkx, blx, bux)
-
-    if len(m.binaries.labels) + len(m.integers.labels) > 0:
-        idx = [i for (i, v) in enumerate(M.vtypes) if v in ["B", "I"]]
-        task.putvartypelist(idx, [mosek.variabletype.type_int] * len(idx))
-        if len(m.binaries.labels) > 0:
-            bidx = [i for (i, v) in enumerate(M.vtypes) if v == "B"]
-            task.putvarboundlistconst(bidx, mosek.boundkey.ra, 0.0, 1.0)
-
-    ## Constraints
-
-    if len(m.constraints) > 0:
-        if set_names:
-            names = print_constraints(M.clabels)
-            for i, n in enumerate(names):
-                task.putconname(i, n)
-        bkc = [
-            (
-                (mosek.boundkey.up if b < np.inf else mosek.boundkey.fr)
-                if s == "<"
-                else (
-                    (mosek.boundkey.lo if b > -np.inf else mosek.boundkey.up)
-                    if s == ">"
-                    else mosek.boundkey.fx
-                )
-            )
-            for s, b in zip(M.sense, M.b)
-        ]
-        blc = [b if b > -np.inf else 0.0 for b in M.b]
-        buc = [b if b < np.inf else 0.0 for b in M.b]
-        # blc = M.b
-        # buc = M.b
-        if M.A is not None:
-            A = M.A.tocsr()
-            task.putarowslice(
-                0, m.ncons, A.indptr[:-1], A.indptr[1:], A.indices, A.data
-            )
-            task.putconboundslice(0, m.ncons, bkc, blc, buc)
-
-    ## Objective
-    if M.Q is not None:
-        Q = (0.5 * tril(M.Q + M.Q.transpose())).tocoo()
-        task.putqobj(Q.row, Q.col, Q.data)
-    task.putclist(list(np.arange(m.nvars)), M.c)
-
-    if m.objective.sense == "max":
-        task.putobjsense(mosek.objsense.maximize)
-    else:
-        task.putobjsense(mosek.objsense.minimize)
-    return task
+    return solvers.Mosek._build_solver_model(
+        m,
+        task,
+        explicit_coordinate_names=explicit_coordinate_names,
+        set_names=set_names,
+    )
 
 
 def to_gurobipy(
@@ -768,84 +659,15 @@ def to_gurobipy(
     explicit_coordinate_names: bool = False,
     set_names: bool = True,
 ) -> Any:
-    """
-    Export the model to gurobipy.
-
-    This function does not write the model to intermediate files but directly
-    passes it to gurobipy. Note that for large models this is not
-    computationally efficient.
-
-    Parameters
-    ----------
-    m : linopy.Model
-    env : gurobipy.Env
-    explicit_coordinate_names : bool, optional
-        Whether to use explicit coordinate names. Default is False.
-    set_names : bool, optional
-        Whether to set variable and constraint names. Default is True.
-        Setting to False can significantly speed up model export.
-
-    Returns
-    -------
-    model : gurobipy.Model
-    """
-    import gurobipy
-
-    m.constraints.sanitize_missings()
-    model = gurobipy.Model(env=env)
-
-    M = m.matrices
-
-    kwargs = {}
-    if set_names:
-        print_variables, print_constraints = get_printers_scalar(
-            m, explicit_coordinate_names=explicit_coordinate_names
-        )
-        kwargs["name"] = print_variables(M.vlabels)
-    if (
-        len(m.binaries.labels)
-        + len(m.integers.labels)
-        + len(list(m.variables.semi_continuous))
-    ):
-        kwargs["vtype"] = M.vtypes
-    x = model.addMVar(M.vlabels.shape, M.lb, M.ub, **kwargs)
-
-    if m.is_quadratic:
-        model.setObjective(0.5 * x.T @ M.Q @ x + M.c @ x)  # type: ignore
-    else:
-        model.setObjective(M.c @ x)
-
-    if m.objective.sense == "max":
-        model.ModelSense = -1
-
-    if len(m.constraints):
-        c = model.addMConstr(M.A, x, M.sense, M.b)  # type: ignore
-        if set_names:
-            names = print_constraints(M.clabels)
-            c.setAttr("ConstrName", names)
-
-    if m.variables.sos:
-        for var_name in m.variables.sos:
-            var = m.variables.sos[var_name]
-            sos_type: int = var.attrs[SOS_TYPE_ATTR]  # type: ignore[assignment]
-            sos_dim: str = var.attrs[SOS_DIM_ATTR]  # type: ignore[assignment]
-
-            def add_sos(s: xr.DataArray, sos_type: int, sos_dim: str) -> None:
-                s = s.squeeze()
-                indices = s.values.flatten().tolist()
-                weights = s.coords[sos_dim].values.tolist()
-                model.addSOS(sos_type, x[indices].tolist(), weights)
-
-            others = [dim for dim in var.labels.dims if dim != sos_dim]
-            if not others:
-                add_sos(var.labels, sos_type, sos_dim)
-            else:
-                stacked = var.labels.stack(_sos_group=others)
-                for _, s in stacked.groupby("_sos_group"):
-                    add_sos(s.unstack("_sos_group"), sos_type, sos_dim)
-
-    model.update()
-    return model
+    """Build the gurobipy.Model for `m`."""
+    solver = solvers.Gurobi.from_model(
+        m,
+        io_api="direct",
+        explicit_coordinate_names=explicit_coordinate_names,
+        set_names=set_names,
+        env=env,
+    )
+    return solver.solver_model
 
 
 def to_highspy(
@@ -853,166 +675,20 @@ def to_highspy(
     explicit_coordinate_names: bool = False,
     set_names: bool = True,
 ) -> Highs:
-    """
-    Export the model to highspy.
-
-    This function does not write the model to intermediate files but directly
-    passes it to highspy.
-
-    Parameters
-    ----------
-    m : linopy.Model
-    explicit_coordinate_names : bool, optional
-        Whether to use explicit coordinate names. Default is False.
-    set_names : bool, optional
-        Whether to set variable and constraint names. Default is True.
-        Setting to False can significantly speed up model export.
-
-    Returns
-    -------
-    model : highspy.Highs
-    """
-    if m.variables.sos:
-        raise NotImplementedError(
-            "SOS constraints are not supported by the HiGHS direct API. "
-            "Use io_api='lp' instead."
-        )
-
-    import highspy
-
-    M = m.matrices
-    h = highspy.Highs()
-    h.addVars(len(M.vlabels), M.lb, M.ub)
-    if len(m.binaries) + len(m.integers) + len(list(m.variables.semi_continuous)):
-        vtypes = M.vtypes
-        # Map linopy vtypes to HiGHS integrality values:
-        # 0 = continuous, 1 = integer, 2 = semi-continuous
-        integrality_map = {"C": 0, "B": 1, "I": 1, "S": 2}
-        int_mask = (vtypes == "B") | (vtypes == "I") | (vtypes == "S")
-        labels = np.arange(len(vtypes))[int_mask]
-        integrality = np.array(
-            [integrality_map[v] for v in vtypes[int_mask]], dtype=np.int32
-        )
-        h.changeColsIntegrality(len(labels), labels, integrality)
-        if len(m.binaries):
-            labels = np.arange(len(vtypes))[vtypes == "B"]
-            n = len(labels)
-            h.changeColsBounds(n, labels, zeros_like(labels), ones_like(labels))
-
-    # linear objective
-    c = M.c
-    h.changeColsCost(len(c), np.arange(len(c), dtype=np.int32), c)
-
-    # linear constraints
-    A = M.A
-    if A is not None:
-        A = A.tocsr()
-        num_cons = A.shape[0]
-        lower = np.where(M.sense != "<", M.b, -np.inf)
-        upper = np.where(M.sense != ">", M.b, np.inf)
-        h.addRows(num_cons, lower, upper, A.nnz, A.indptr, A.indices, A.data)
-
-    if set_names:
-        print_variables, print_constraints = get_printers_scalar(
-            m, explicit_coordinate_names=explicit_coordinate_names
-        )
-        lp = h.getLp()
-        lp.col_names_ = print_variables(M.vlabels)
-        if len(M.clabels):
-            lp.row_names_ = print_constraints(M.clabels)
-        h.passModel(lp)
-
-    # quadrative objective
-    Q = M.Q
-    if Q is not None:
-        Q = triu(Q)
-        Q = Q.tocsr()
-        num_vars = Q.shape[0]
-        h.passHessian(num_vars, Q.nnz, 1, Q.indptr, Q.indices, Q.data)
-
-    # change objective sense
-    if m.objective.sense == "max":
-        h.changeObjectiveSense(highspy.ObjSense.kMaximize)
-
-    return h
-
-
-def to_cupdlpx(
-    m: Model,
-    explicit_coordinate_names: bool = False,
-    set_names: bool = True,
-) -> cupdlpxModel:
-    """
-    Export the model to cupdlpx.
-
-    This function does not write the model to intermediate files but directly
-    passes it to cupdlpx.
-
-    cuPDLPx does not support named variables and constraints, so the
-    `explicit_coordinate_names` and `set_names` parameters are ignored.
-
-    Parameters
-    ----------
-    m : linopy.Model
-    explicit_coordinate_names : bool, optional
-        Ignored. cuPDLPx does not support named variables/constraints.
-    set_names : bool, optional
-        Ignored. cuPDLPx does not support named variables/constraints.
-
-    Returns
-    -------
-    model : cupdlpx.Model
-    """
-    if m.variables.semi_continuous:
-        raise NotImplementedError(
-            "Semi-continuous variables are not supported by cuPDLPx. "
-            "Use a solver that supports them (gurobi, cplex, highs)."
-        )
-
-    import cupdlpx
-
-    if explicit_coordinate_names:
-        warnings.warn(
-            "cuPDLPx does not support named variables/constraints. "
-            "The explicit_coordinate_names parameter is ignored.",
-            UserWarning,
-            stacklevel=2,
-        )
-
-    # build model using canonical form matrices and vectors
-    # see https://github.com/MIT-Lu-Lab/cuPDLPx/tree/main/python#modeling
-    M = m.matrices
-    if M.A is None:
-        msg = "Model has no constraints, cannot export to cuPDLPx."
-        raise ValueError(msg)
-    A = M.A.tocsr()  # cuPDLPx only supports CSR sparse matrix format
-    # linopy stores constraints as Ax ?= b and keeps track of inequality
-    # sense in M.sense. Convert to separate lower and upper bound vectors.
-    l = np.where(
-        np.logical_or(np.equal(M.sense, ">"), np.equal(M.sense, "=")),
-        M.b,
-        -np.inf,
-    )
-    u = np.where(
-        np.logical_or(np.equal(M.sense, "<"), np.equal(M.sense, "=")),
-        M.b,
-        np.inf,
-    )
-
-    cu_model = cupdlpx.Model(
-        objective_vector=M.c,
-        constraint_matrix=A,
-        constraint_lower_bound=l,
-        constraint_upper_bound=u,
-        variable_lower_bound=M.lb,
-        variable_upper_bound=M.ub,
+    """Build the highspy.Highs instance for `m`."""
+    solver = solvers.Highs.from_model(
+        m,
+        io_api="direct",
+        explicit_coordinate_names=explicit_coordinate_names,
+        set_names=set_names,
     )
+    return solver.solver_model
 
-    # change objective sense
-    if m.objective.sense == "max":
-        cu_model.ModelSense = cupdlpx.PDLP.MAXIMIZE
 
-    return cu_model
+def to_cupdlpx(m: Model) -> cupdlpxModel:
+    """Build the cupdlpx.Model for `m`."""
+    solver = solvers.cuPDLPx.from_model(m, io_api="direct")
+    return solver.solver_model
 
 
 def to_block_files(m: Model, fn: Path) -> None:
diff --git a/linopy/matrices.py b/linopy/matrices.py
index 1fb59344..e694a720 100644
--- a/linopy/matrices.py
+++ b/linopy/matrices.py
@@ -81,18 +81,16 @@ def _build_cons(self) -> None:
 
         label_index = m.variables.label_index
         csrs = []
-        clabels_list = []
         b_list = []
         sense_list = []
         for c in m.constraints.data.values():
-            csr, con_labels, b, sense = c.to_matrix_with_rhs(label_index)
+            csr, _, b, sense = c.to_matrix_with_rhs(label_index)
             csrs.append(csr)
-            clabels_list.append(con_labels)
             b_list.append(b)
             sense_list.append(sense)
 
         self.A = cast(scipy.sparse.csr_array, scipy.sparse.vstack(csrs, format="csr"))
-        self.clabels = np.concatenate(clabels_list)
+        self.clabels = m.constraints.label_index.clabels
         self.b = np.concatenate(b_list) if b_list else np.array([])
         self.sense = (
             np.concatenate(sense_list) if sense_list else np.array([], dtype=object)
diff --git a/linopy/model.py b/linopy/model.py
index 21e4e29c..665cad84 100644
--- a/linopy/model.py
+++ b/linopy/model.py
@@ -20,7 +20,7 @@
 import pandas as pd
 import xarray as xr
 from deprecation import deprecated
-from numpy import inf, nan, ndarray
+from numpy import inf, ndarray
 from pandas.core.frame import DataFrame
 from pandas.core.series import Series
 from xarray import DataArray, Dataset
@@ -32,11 +32,8 @@
     assign_multiindex_safe,
     best_int,
     broadcast_mask,
-    lookup_vals,
     maybe_replace_signs,
     replace_by_map,
-    series_to_lookup_array,
-    set_int_index,
     to_path,
 )
 from linopy.constants import (
@@ -48,6 +45,7 @@
     SOS_TYPE_ATTR,
     TERM_DIM,
     ModelStatus,
+    Result,
     TerminationCondition,
 )
 from linopy.constraints import (
@@ -85,9 +83,9 @@
     from linopy.remote import OetcHandler
 except ImportError:
     OetcHandler = None  # type: ignore
-from linopy.solver_capabilities import SolverFeature, solver_supports
 from linopy.solvers import (
     IO_APIS,
+    SolverFeature,
     available_solvers,
 )
 from linopy.sos_reformulation import (
@@ -193,8 +191,7 @@ class Model:
     the optimization process.
     """
 
-    solver_model: Any
-    solver_name: str
+    _solver: solvers.Solver | None
     _variables: Variables
     _constraints: Constraints
     _objective: Objective
@@ -241,8 +238,7 @@ class Model:
         "_solver_dir",
         "_relaxed_registry",
         "_piecewise_formulations",
-        "solver_model",
-        "solver_name",
+        "_solver",
         "__weakref__",
     )
 
@@ -312,6 +308,37 @@ def __init__(
         self._solver_dir: Path = Path(
             gettempdir() if solver_dir is None else solver_dir
         )
+        self._solver: solvers.Solver | None = None
+
+    @property
+    def solver(self) -> solvers.Solver | None:
+        return self._solver
+
+    @solver.setter
+    def solver(self, value: solvers.Solver | None) -> None:
+        if self._solver is not None and self._solver is not value:
+            self._solver.close()
+        self._solver = value
+
+    @property
+    def solver_model(self) -> Any:
+        return self.solver.solver_model if self.solver is not None else None
+
+    @solver_model.setter
+    def solver_model(self, value: Any) -> None:
+        if value is not None:
+            raise AttributeError("solver state is managed via model.solver")
+        self.solver = None
+
+    @property
+    def solver_name(self) -> str | None:
+        return self.solver.solver_name.value if self.solver is not None else None
+
+    @solver_name.setter
+    def solver_name(self, value: str | None) -> None:
+        if value is not None:
+            raise AttributeError("solver state is managed via model.solver")
+        self.solver = None
 
     @property
     def matrices(self) -> MatrixAccessor:
@@ -1648,11 +1675,13 @@ def solve(
             )
             logger.info(f"Solver options:\n{options_string}")
 
+        solver_class = getattr(solvers, solvers.SolverName(solver_name).name)
+
         if problem_fn is None:
             problem_fn = self.get_problem_file(io_api=io_api)
         if solution_fn is None:
             if (
-                solver_supports(solver_name, SolverFeature.SOLUTION_FILE_NOT_NEEDED)
+                solver_class.supports(SolverFeature.SOLUTION_FILE_NOT_NEEDED)
                 and not keep_files
             ):
                 # these (solver, keep_files=False) combos do not need a solution file
@@ -1666,8 +1695,8 @@ def solve(
         if sanitize_infinities:
             self.constraints.sanitize_infinities()
 
-        if self.is_quadratic and not solver_supports(
-            solver_name, SolverFeature.QUADRATIC_OBJECTIVE
+        if self.is_quadratic and not solver_class.supports(
+            SolverFeature.QUADRATIC_OBJECTIVE
         ):
             raise ValueError(
                 f"Solver {solver_name} does not support quadratic problems."
@@ -1681,7 +1710,7 @@ def solve(
 
         sos_reform_result = None
         if self.variables.sos:
-            supports_sos = solver_supports(solver_name, SolverFeature.SOS_CONSTRAINTS)
+            supports_sos = solver_class.supports(SolverFeature.SOS_CONSTRAINTS)
             if reformulate_sos in (True, "auto") and not supports_sos:
                 logger.info(f"Reformulating SOS constraints for solver {solver_name}")
                 sos_reform_result = reformulate_sos_constraints(self)
@@ -1697,107 +1726,92 @@ def solve(
                 )
 
         if self.variables.semi_continuous:
-            if not solver_supports(
-                solver_name, SolverFeature.SEMI_CONTINUOUS_VARIABLES
-            ):
+            if not solver_class.supports(SolverFeature.SEMI_CONTINUOUS_VARIABLES):
                 raise ValueError(
                     f"Solver {solver_name} does not support semi-continuous variables. "
                     "Use a solver that supports them (gurobi, cplex, highs)."
                 )
 
         try:
-            solver_class = getattr(solvers, f"{solvers.SolverName(solver_name).name}")
-            # initialize the solver as object of solver subclass <solver_class>
-            solver = solver_class(
-                **solver_options,
-            )
+            self.solver = None  # closes any previous solver
             if io_api == "direct":
                 if set_names is None:
                     set_names = self.set_names_in_solver_io
-                # no problem file written and direct model is set for solver
-                result = solver.solve_problem_from_model(
-                    model=self,
-                    solution_fn=to_path(solution_fn),
-                    log_fn=to_path(log_fn),
-                    warmstart_fn=to_path(warmstart_fn),
-                    basis_fn=to_path(basis_fn),
-                    env=env,
-                    explicit_coordinate_names=explicit_coordinate_names,
-                    set_names=set_names,
-                )
+                build_kwargs: dict[str, Any] = {
+                    "explicit_coordinate_names": explicit_coordinate_names,
+                    "set_names": set_names,
+                    "log_fn": to_path(log_fn),
+                }
+                if env is not None:
+                    build_kwargs["env"] = env
             else:
-                if (
-                    not solver_supports(solver_name, SolverFeature.LP_FILE_NAMES)
-                    and explicit_coordinate_names
-                ):
-                    logger.warning(
-                        f"{solver_name} does not support writing names to lp files, disabling it."
-                    )
-                    explicit_coordinate_names = False
-                problem_fn = self.to_file(
-                    to_path(problem_fn),
-                    io_api=io_api,
-                    explicit_coordinate_names=explicit_coordinate_names,
-                    slice_size=slice_size,
-                    progress=progress,
-                )
-                result = solver.solve_problem_from_file(
-                    problem_fn=to_path(problem_fn),
-                    solution_fn=to_path(solution_fn),
-                    log_fn=to_path(log_fn),
-                    warmstart_fn=to_path(warmstart_fn),
-                    basis_fn=to_path(basis_fn),
-                    env=env,
-                )
-
+                build_kwargs = {
+                    "explicit_coordinate_names": explicit_coordinate_names,
+                    "slice_size": slice_size,
+                    "progress": progress,
+                    "problem_fn": to_path(problem_fn),
+                }
+            self.solver = solver = solvers.Solver.from_name(
+                solver_name,
+                model=self,
+                io_api=io_api,
+                options=solver_options,
+                **build_kwargs,
+            )
+            if io_api != "direct":
+                problem_fn = solver._problem_fn
+            result = solver.solve(
+                solution_fn=to_path(solution_fn),
+                log_fn=to_path(log_fn),
+                warmstart_fn=to_path(warmstart_fn),
+                basis_fn=to_path(basis_fn),
+                env=env,
+            )
         finally:
             for fn in (problem_fn, solution_fn):
                 if fn is not None and (os.path.exists(fn) and not keep_files):
                     os.remove(fn)
 
         try:
-            result.info()
+            return self.apply_result(result)
+        finally:
+            if sos_reform_result is not None:
+                undo_sos_reformulation(self, sos_reform_result)
 
-            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,
-                )
+    def apply_result(self, result: Result) -> tuple[str, str]:
+        result.info()
 
-            # 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
+        if result.solution is not None:
+            self.objective._value = result.solution.objective
 
-            sol_arr = series_to_lookup_array(sol)
+        status_value = result.status.status.value
+        termination_condition = result.status.termination_condition.value
+        self.status = status_value
+        self.termination_condition = termination_condition
 
-            for _, var in self.variables.items():
-                vals = lookup_vals(sol_arr, np.ravel(var.labels))
-                var.solution = xr.DataArray(vals.reshape(var.labels.shape), var.coords)
+        if not result.status.is_ok:
+            return status_value, termination_condition
 
-            if not result.solution.dual.empty:
-                dual = result.solution.dual.copy()
-                dual = set_int_index(dual)
-                dual.loc[-1] = nan
+        if result.solution is None or len(result.solution.primal) == 0:
+            return status_value, termination_condition
 
-                dual_arr = series_to_lookup_array(dual)
+        primal = result.solution.primal
+        for _, var in self.variables.items():
+            start, end = var.range
+            var.solution = xr.DataArray(
+                primal[start:end].reshape(var.shape), var.coords
+            )
 
-                for _, con in self.constraints.items():
-                    vals = lookup_vals(dual_arr, np.ravel(con.labels))
-                    con.dual = xr.DataArray(
-                        vals.reshape(con.labels.shape), con.labels.coords
-                    )
+        if len(result.solution.dual):
+            dual = result.solution.dual
+            for _, con in self.constraints.items():
+                start, end = con.range
+                coords = {dim: con.coords[dim] for dim in con.coord_dims}
+                con.dual = xr.DataArray(
+                    dual[start:end].reshape(con.shape), coords, dims=con.coord_dims
+                )
 
-            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)
+        return status_value, termination_condition
 
     def _mock_solve(
         self,
@@ -1807,6 +1821,7 @@ def _mock_solve(
         solver_name = "mock"
 
         logger.info(f" Solve problem using {solver_name.title()} solver")
+        self.solver = None
         # reset result
         self.reset_solution()
 
@@ -1819,8 +1834,6 @@ def _mock_solve(
         self.objective._value = 0.0
         self.status = "ok"
         self.termination_condition = TerminationCondition.optimal.value
-        self.solver_model = None
-        self.solver_name = solver_name
 
         for name, var in self.variables.items():
             var.solution = xr.DataArray(0.0, var.coords)
@@ -1843,7 +1856,7 @@ def compute_infeasibilities(self) -> list[int]:
         labels : list[int]
             Labels of the infeasible constraints.
         """
-        solver_model = getattr(self, "solver_model", None)
+        solver_model = self.solver_model
 
         # Check for Gurobi
         if "gurobi" in available_solvers:
@@ -1872,8 +1885,10 @@ def compute_infeasibilities(self) -> list[int]:
         # If we get here, either the solver doesn't support IIS or no solver model is available
         if solver_model is None:
             # Check if this is a supported solver without a stored model
-            solver_name = getattr(self, "solver_name", "unknown")
-            if solver_supports(solver_name, SolverFeature.IIS_COMPUTATION):
+            solver_name = self.solver_name or "unknown"
+            if self.solver is not None and self.solver.supports(
+                SolverFeature.IIS_COMPUTATION
+            ):
                 raise ValueError(
                     "No solver model available. The model must be solved first with "
                     "a solver that supports IIS computation and the result must be infeasible."
@@ -1930,7 +1945,7 @@ def _compute_infeasibilities_xpress(self, solver_model: Any) -> list[int]:
         if solver_model.attributes.numiis == 0:
             return []
 
-        clabels = self.matrices.clabels
+        clabels = self.constraints.label_index.clabels
         constraint_position_map = {}
         for position, constraint_obj in enumerate(solver_model.getConstraint()):
             if 0 <= position < len(clabels):
diff --git a/linopy/solver_capabilities.py b/linopy/solver_capabilities.py
index f9c6aba4..0e748082 100644
--- a/linopy/solver_capabilities.py
+++ b/linopy/solver_capabilities.py
@@ -1,306 +1,100 @@
 """
-Linopy module for solver capability tracking.
+Back-compat shim for legacy solver-capability imports.
 
-This module provides a centralized registry of solver capabilities,
-replacing scattered hardcoded checks throughout the codebase.
+Capability data is declared on each `Solver` subclass in `linopy.solvers`.
+Prefer `Solver.features` / `Solver.supports()` over the helpers in this module.
 """
 
 from __future__ import annotations
 
+from collections.abc import Iterator, Mapping, Sequence
 from dataclasses import dataclass
-from enum import Enum, auto
-from importlib.metadata import PackageNotFoundError
-from importlib.metadata import version as package_version
+from enum import Enum
 from typing import TYPE_CHECKING
 
-from packaging.specifiers import SpecifierSet
-
 if TYPE_CHECKING:
-    from collections.abc import Sequence
-
-
-def _xpress_supports_gpu() -> bool:
-    """Check if installed xpress version supports GPU acceleration (>=9.8.0)."""
-    try:
-        return package_version("xpress") in SpecifierSet(">=9.8.0")
-    except PackageNotFoundError:
-        return False
-
-
-class SolverFeature(Enum):
-    """Enumeration of all solver capabilities tracked by linopy."""
+    from linopy.solvers import Solver, SolverFeature
 
-    # Model feature support
-    INTEGER_VARIABLES = auto()  # Support for integer variables
+__all__ = (
+    "SOLVER_REGISTRY",
+    "SolverFeature",
+    "SolverInfo",
+    "get_available_solvers_with_feature",
+    "get_solvers_with_feature",
+    "solver_supports",
+)
 
-    # Objective function support
-    QUADRATIC_OBJECTIVE = auto()
 
-    # I/O capabilities
-    DIRECT_API = auto()  # Solve directly from Model without writing files
-    LP_FILE_NAMES = auto()  # Support for named variables/constraints in LP files
-    READ_MODEL_FROM_FILE = auto()  # Ability to read models from file
-    SOLUTION_FILE_NOT_NEEDED = auto()  # Solver doesn't need a solution file
+def __getattr__(name: str) -> object:
+    if name == "SolverFeature":
+        from linopy import solvers as _solvers_mod
 
-    # Advanced features
-    GPU_ACCELERATION = auto()  # GPU-accelerated solving
-    IIS_COMPUTATION = auto()  # Irreducible Infeasible Set computation
-
-    # Special constraint types
-    SOS_CONSTRAINTS = auto()  # Special Ordered Sets (SOS1/SOS2) constraints
-
-    # Special variable types
-    SEMI_CONTINUOUS_VARIABLES = auto()  # Semi-continuous variable support
-
-    # Solver-specific
-    SOLVER_ATTRIBUTE_ACCESS = auto()  # Direct access to solver variable attributes
+        return _solvers_mod.SolverFeature
+    raise AttributeError(name)
 
 
 @dataclass(frozen=True)
 class SolverInfo:
-    """Information about a solver's capabilities."""
+    """Legacy view of a solver's capabilities. Prefer Solver.features / Solver.supports()."""
 
     name: str
-    features: frozenset[SolverFeature]
+    features: frozenset[Enum]
     display_name: str = ""
 
     def __post_init__(self) -> None:
         if not self.display_name:
             object.__setattr__(self, "display_name", self.name.upper())
 
-    def supports(self, feature: SolverFeature) -> bool:
-        """Check if this solver supports a given feature."""
+    def supports(self, feature: Enum) -> bool:
         return feature in self.features
 
 
-# Define all solver capabilities
-SOLVER_REGISTRY: dict[str, SolverInfo] = {
-    "gurobi": SolverInfo(
-        name="gurobi",
-        display_name="Gurobi",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.DIRECT_API,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-                SolverFeature.IIS_COMPUTATION,
-                SolverFeature.SOS_CONSTRAINTS,
-                SolverFeature.SEMI_CONTINUOUS_VARIABLES,
-                SolverFeature.SOLVER_ATTRIBUTE_ACCESS,
-            }
-        ),
-    ),
-    "highs": SolverInfo(
-        name="highs",
-        display_name="HiGHS",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.DIRECT_API,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-                SolverFeature.SEMI_CONTINUOUS_VARIABLES,
-            }
-        ),
-    ),
-    "glpk": SolverInfo(
-        name="glpk",
-        display_name="GLPK",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-            }
-        ),  # No LP_FILE_NAMES support
-    ),
-    "cbc": SolverInfo(
-        name="cbc",
-        display_name="CBC",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-            }
-        ),  # No LP_FILE_NAMES support
-    ),
-    "cplex": SolverInfo(
-        name="cplex",
-        display_name="CPLEX",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOS_CONSTRAINTS,
-                SolverFeature.SEMI_CONTINUOUS_VARIABLES,
-            }
-        ),
-    ),
-    "xpress": SolverInfo(
-        name="xpress",
-        display_name="FICO Xpress",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-                SolverFeature.GPU_ACCELERATION,
-                SolverFeature.IIS_COMPUTATION,
-            }
-            if _xpress_supports_gpu()
-            else {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-                SolverFeature.IIS_COMPUTATION,
-            }
-        ),
-    ),
-    "knitro": SolverInfo(
-        name="knitro",
-        display_name="Artelys Knitro",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-            }
-        ),
-    ),
-    "scip": SolverInfo(
-        name="scip",
-        display_name="SCIP",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-            }
-        ),
-    ),
-    "mosek": SolverInfo(
-        name="mosek",
-        display_name="MOSEK",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.DIRECT_API,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-            }
-        ),
-    ),
-    "copt": SolverInfo(
-        name="copt",
-        display_name="COPT",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-            }
-        ),
-    ),
-    "mindopt": SolverInfo(
-        name="mindopt",
-        display_name="MindOpt",
-        features=frozenset(
-            {
-                SolverFeature.INTEGER_VARIABLES,
-                SolverFeature.QUADRATIC_OBJECTIVE,
-                SolverFeature.LP_FILE_NAMES,
-                SolverFeature.READ_MODEL_FROM_FILE,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-            }
-        ),
-    ),
-    "cupdlpx": SolverInfo(
-        name="cupdlpx",
-        display_name="cuPDLPx",
-        features=frozenset(
-            {
-                SolverFeature.DIRECT_API,
-                SolverFeature.GPU_ACCELERATION,
-                SolverFeature.SOLUTION_FILE_NOT_NEEDED,
-            }
-        ),
-    ),
-}
+def _solver_class(name: str) -> type[Solver] | None:
+    from linopy import solvers as _solvers_mod
 
+    try:
+        return getattr(_solvers_mod, _solvers_mod.SolverName(name).name, None)
+    except ValueError:
+        return None
 
-def solver_supports(solver_name: str, feature: SolverFeature) -> bool:
-    """
-    Check if a solver supports a given feature.
-
-    Parameters
-    ----------
-    solver_name : str
-        Name of the solver (e.g., "gurobi", "highs")
-    feature : SolverFeature
-        The feature to check for
 
-    Returns
-    -------
-    bool
-        True if the solver supports the feature, False otherwise.
-        Returns False for unknown solvers.
-    """
-    if solver_name not in SOLVER_REGISTRY:
-        return False
-    return SOLVER_REGISTRY[solver_name].supports(feature)
+def solver_supports(solver_name: str, feature: SolverFeature) -> bool:
+    cls = _solver_class(solver_name)
+    return cls is not None and cls.supports(feature)
 
 
 def get_solvers_with_feature(feature: SolverFeature) -> list[str]:
-    """
-    Get all solvers that support a given feature.
-
-    Parameters
-    ----------
-    feature : SolverFeature
-        The feature to filter by
+    from linopy.solvers import SolverName
 
-    Returns
-    -------
-    list[str]
-        List of solver names supporting the feature
-    """
-    return [name for name, info in SOLVER_REGISTRY.items() if info.supports(feature)]
+    return [n.value for n in SolverName if solver_supports(n.value, feature)]
 
 
 def get_available_solvers_with_feature(
     feature: SolverFeature, available_solvers: Sequence[str]
 ) -> list[str]:
-    """
-    Get installed solvers that support a given feature.
+    return [s for s in get_solvers_with_feature(feature) if s in available_solvers]
 
-    Parameters
-    ----------
-    feature : SolverFeature
-        The feature to filter by
-    available_solvers : Sequence[str]
-        List of currently available/installed solvers
 
-    Returns
-    -------
-    list[str]
-        List of installed solver names supporting the feature
-    """
-    return [s for s in get_solvers_with_feature(feature) if s in available_solvers]
+class _LazyRegistry(Mapping[str, SolverInfo]):
+    def __getitem__(self, key: str) -> SolverInfo:
+        cls = _solver_class(key)
+        if cls is None:
+            raise KeyError(key)
+        return SolverInfo(
+            name=key,
+            features=cls.supported_features(),
+            display_name=cls.display_name,
+        )
+
+    def __iter__(self) -> Iterator[str]:
+        from linopy.solvers import SolverName
+
+        return (n.value for n in SolverName)
+
+    def __len__(self) -> int:
+        from linopy.solvers import SolverName
+
+        return len(SolverName)
+
+
+SOLVER_REGISTRY: Mapping[str, SolverInfo] = _LazyRegistry()
diff --git a/linopy/solvers.py b/linopy/solvers.py
index 86c312e4..e76691a4 100644
--- a/linopy/solvers.py
+++ b/linopy/solvers.py
@@ -15,42 +15,130 @@
 import sys
 import threading
 import warnings
-from abc import ABC, abstractmethod
+from abc import ABC
 from collections import namedtuple
 from collections.abc import Callable, Generator
+from dataclasses import dataclass, field
+from enum import Enum, auto
+from importlib.metadata import PackageNotFoundError
+from importlib.metadata import version as package_version
 from pathlib import Path
-from typing import TYPE_CHECKING, Any, Generic, TypeVar
+from typing import TYPE_CHECKING, Any, ClassVar, Generic, TypeVar
 
 import numpy as np
 import pandas as pd
+import xarray as xr
+from packaging.specifiers import SpecifierSet
 from packaging.version import parse as parse_version
+from scipy.sparse import tril, triu
 
 import linopy.io
+from linopy.common import count_initial_letters, values_to_lookup_array
 from linopy.constants import (
+    SOS_DIM_ATTR,
+    SOS_TYPE_ATTR,
     Result,
     Solution,
+    SolverReport,
     SolverStatus,
     Status,
     TerminationCondition,
 )
-from linopy.solver_capabilities import (
-    SolverFeature,
-    get_solvers_with_feature,
-)
+
+
+def _parse_int_label(name: str) -> int:
+    """Strip leading non-digits and parse the integer label."""
+    s = str(name)
+    cutoff = count_initial_letters(s)
+    try:
+        return int(s[cutoff:])
+    except ValueError:
+        return int(re.sub(r".*#", "", s))
+
+
+def _names_to_labels(names: Any) -> np.ndarray:
+    """Vectorised conversion of solver-provided names to integer labels."""
+    if len(names) == 0:
+        return np.array([], dtype=np.int64)
+    index = pd.Index(names)
+    if pd.api.types.is_integer_dtype(index.dtype):
+        return index.to_numpy(dtype=np.int64)
+    string_index = index.astype(str)
+    cutoff = count_initial_letters(str(string_index[0]))
+    try:
+        return string_index.str[cutoff:].astype(np.int64).to_numpy(dtype=np.int64)
+    except (TypeError, ValueError):
+        try:
+            return (
+                string_index.str.replace(r".*#", "", regex=True)
+                .astype(np.int64)
+                .to_numpy(dtype=np.int64)
+            )
+        except (TypeError, ValueError):
+            return np.fromiter(
+                (_parse_int_label(n) for n in names), dtype=np.int64, count=len(names)
+            )
+
+
+def _solution_from_names(values: np.ndarray, names: Any, size: int) -> np.ndarray:
+    """
+    Build a label-indexed dense solution array of length ``size`` from
+    solver-side names. Used by paths where the solver may iterate in arbitrary
+    order or drop unused entities (file-based LP solvers, the ``from_file``
+    paths of Highs/Gurobi).
+    """
+    if not size:
+        return np.array([], dtype=float)
+    return values_to_lookup_array(
+        np.asarray(values, dtype=float), _names_to_labels(names), size=size
+    )
+
+
+def _solution_from_labels(
+    values: np.ndarray, labels: np.ndarray | None, size: int
+) -> np.ndarray:
+    """Scatter solver-side values into a label-indexed dense array of length ``size``."""
+    if not size:
+        return np.array([], dtype=float)
+    assert labels is not None
+    return values_to_lookup_array(np.asarray(values, dtype=float), labels, size=size)
+
+
+class SolverFeature(Enum):
+    """Enumeration of all solver capabilities tracked by linopy."""
+
+    INTEGER_VARIABLES = auto()
+    QUADRATIC_OBJECTIVE = auto()
+    DIRECT_API = auto()
+    LP_FILE_NAMES = auto()
+    READ_MODEL_FROM_FILE = auto()
+    SOLUTION_FILE_NOT_NEEDED = auto()
+    GPU_ACCELERATION = auto()
+    IIS_COMPUTATION = auto()
+    SOS_CONSTRAINTS = auto()
+    SEMI_CONTINUOUS_VARIABLES = auto()
+    SOLVER_ATTRIBUTE_ACCESS = auto()
+    MIP_DUAL_BOUND_REPORT = auto()
+
+
+def _installed_version_in(pkg: str, spec: str) -> bool:
+    """Check whether the installed version of `pkg` satisfies `spec`."""
+    try:
+        return package_version(pkg) in SpecifierSet(spec)
+    except PackageNotFoundError:
+        return False
+
 
 if TYPE_CHECKING:
+    import cupdlpx
     import gurobipy
+    import highspy
+    import mosek
 
     from linopy.model import Model
 
 EnvType = TypeVar("EnvType")
 
-# Generated from solver_capabilities registry for backward compatibility
-QUADRATIC_SOLVERS = get_solvers_with_feature(SolverFeature.QUADRATIC_OBJECTIVE)
-NO_SOLUTION_FILE_SOLVERS = get_solvers_with_feature(
-    SolverFeature.SOLUTION_FILE_NOT_NEEDED
-)
-
 FILE_IO_APIS = ["lp", "lp-polars", "mps"]
 IO_APIS = FILE_IO_APIS + ["direct"]
 
@@ -220,7 +308,6 @@ class xpress_Namespaces:  # type: ignore[no-redef]
         pass
 
 
-quadratic_solvers = [s for s in QUADRATIC_SOLVERS if s in available_solvers]
 logger = logging.getLogger(__name__)
 
 
@@ -297,84 +384,173 @@ def maybe_adjust_objective_sign(
     return solution
 
 
+@dataclass
 class Solver(ABC, Generic[EnvType]):
     """
     Abstract base class for solving a given linear problem.
 
-    All relevant functions are passed on to the specific solver subclasses.
-    Subclasses must implement the `solve_problem_from_model()` and
-    `solve_problem_from_file()` methods.
+    Subclasses provide ``_build_direct`` / ``_run_direct`` (when supporting the
+    direct API) and ``_run_file`` (when supporting LP/MPS files). Construction
+    goes via :meth:`Solver.from_name` or :meth:`Solver.from_model`.
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        self.solver_options = solver_options
-
-        # Check for the solver to be initialized whether the package is installed or not.
+    model: Model | None = None
+    io_api: str | None = None
+    options: dict[str, Any] = field(default_factory=dict)
+
+    # Runtime state — never set via constructor.
+    status: Status | None = field(init=False, default=None, repr=False)
+    solution: Solution | None = field(init=False, default=None, repr=False)
+    report: SolverReport | None = field(init=False, default=None, repr=False)
+    solver_model: Any = field(init=False, default=None, repr=False)
+    sense: str | None = field(init=False, default=None, repr=False)
+    env: Any = field(init=False, default=None, repr=False)
+    _env_stack: contextlib.ExitStack | None = field(
+        init=False, default=None, repr=False
+    )
+    _vlabels: np.ndarray | None = field(init=False, default=None, repr=False)
+    _clabels: np.ndarray | None = field(init=False, default=None, repr=False)
+    _n_vars: int = field(init=False, default=0, repr=False)
+    _n_cons: int = field(init=False, default=0, repr=False)
+    _problem_fn: Path | None = field(init=False, default=None, repr=False)
+
+    display_name: ClassVar[str] = ""
+    features: ClassVar[frozenset[SolverFeature]] = frozenset()
+    accepted_io_apis: ClassVar[frozenset[str]] = frozenset()
+
+    def __post_init__(self) -> None:
+        if type(self) is Solver:
+            raise TypeError(
+                "Solver is abstract; instantiate a concrete subclass instead."
+            )
         if self.solver_name.value not in available_solvers:
-            msg = f"Solver package for '{self.solver_name.value}' is not installed. Please install first to initialize solver instance."
+            msg = (
+                f"Solver package for '{self.solver_name.value}' is not installed. "
+                "Please install first to initialize solver instance."
+            )
             raise ImportError(msg)
 
-    def safe_get_solution(
-        self, status: Status, func: Callable[[], Solution]
-    ) -> Solution:
+    @property
+    def solver_options(self) -> dict[str, Any]:
+        """Back-compat alias for ``self.options``."""
+        return self.options
+
+    @classmethod
+    def runtime_features(cls) -> frozenset[SolverFeature]:
         """
-        Get solution from function call, if status is unknown still try to run it.
+        Features whose availability depends on the installed solver version
+        or runtime environment. Override in subclasses; the default is empty.
         """
-        if status.is_ok:
-            return func()
-        elif status.status == SolverStatus.unknown:
-            try:
-                logger.warning("Solution status unknown. Trying to parse solution.")
-                sol = func()
-                status.status = SolverStatus.ok
-                logger.warning("Solution parsed successfully.")
-                return sol
-            except Exception as e:
-                logger.error(f"Failed to parse solution: {e}")
-        return Solution()
+        return frozenset()
 
-    @abstractmethod
-    def solve_problem_from_model(
-        self,
+    @classmethod
+    def supported_features(cls) -> frozenset[SolverFeature]:
+        """All features supported by this solver, static plus runtime."""
+        return cls.features | cls.runtime_features()
+
+    @classmethod
+    def supports(cls, feature: SolverFeature) -> bool:
+        """Check if this solver supports a given feature."""
+        return feature in cls.features or feature in cls.runtime_features()
+
+    @staticmethod
+    def from_name(
+        name: str,
         model: Model,
-        solution_fn: Path | None = None,
-        log_fn: Path | None = None,
-        warmstart_fn: Path | None = None,
-        basis_fn: Path | None = None,
-        env: EnvType | None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        """
-        Abstract method to solve a linear problem from a model.
+        io_api: str | None = None,
+        options: dict[str, Any] | None = None,
+        **build_kwargs: Any,
+    ) -> Solver:
+        """Construct and build the solver subclass registered as ``name``."""
+        cls = _solver_class_for(name)
+        if cls is None:
+            raise ValueError(f"unknown solver: {name}")
+        return cls.from_model(
+            model, io_api=io_api, options=options or {}, **build_kwargs
+        )
 
-        Needs to be implemented in the specific solver subclass. Even if the solver
-        does not support solving from a model, this method should be implemented and
-        raise a NotImplementedError.
-        """
-        pass
+    @classmethod
+    def from_model(
+        cls,
+        model: Model,
+        io_api: str | None = None,
+        options: dict[str, Any] | None = None,
+        **build_kwargs: Any,
+    ) -> Solver:
+        """Instantiate and build the solver against ``model``."""
+        instance = cls(model=model, io_api=io_api, options=options or {})
+        instance._build(**build_kwargs)
+        return instance
+
+    def _build(self, **build_kwargs: Any) -> None:
+        """Dispatch to direct or file build based on ``io_api``."""
+        if self.model is None:
+            raise RuntimeError("Solver has no model attached; cannot build.")
+        if self.io_api == "direct":
+            self._build_direct(**build_kwargs)
+        else:
+            self._build_file(**build_kwargs)
 
-    @abstractmethod
-    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: EnvType | None = None,
-    ) -> Result:
-        """
-        Abstract method to solve a linear problem from a problem file.
+    def _build_direct(self, **build_kwargs: Any) -> None:
+        """Build the native solver model from ``self.model``. Override per-solver."""
+        raise NotImplementedError(
+            f"Solver {self.solver_name.value} does not support direct API model export."
+        )
 
-        Needs to be implemented in the specific solver subclass. Even if the solver
-        does not support solving from a file, this method should be implemented and
-        raise a NotImplementedError.
-        """
-        pass
+    def _build_file(self, **build_kwargs: Any) -> None:
+        """Write the LP/MPS file for ``self.model`` and cache its path."""
+        model = self.model
+        assert model is not None
+        io_api = self.io_api
+        if io_api is not None and io_api not in FILE_IO_APIS:
+            raise ValueError(
+                f"Keyword argument `io_api` has to be one of {IO_APIS} or None"
+            )
+        explicit_coordinate_names = build_kwargs.pop("explicit_coordinate_names", False)
+        slice_size = build_kwargs.pop("slice_size", 2_000_000)
+        progress = build_kwargs.pop("progress", None)
+        problem_fn = build_kwargs.pop("problem_fn", None)
+        if problem_fn is None:
+            problem_fn = model.get_problem_file(io_api=io_api)
+        if not self.supports(SolverFeature.LP_FILE_NAMES) and explicit_coordinate_names:
+            logger.warning(
+                f"{self.solver_name.value} does not support writing names to "
+                "lp files, disabling it."
+            )
+            explicit_coordinate_names = False
+        problem_fn = model.to_file(
+            Path(problem_fn) if not isinstance(problem_fn, Path) else problem_fn,
+            io_api=io_api,
+            explicit_coordinate_names=explicit_coordinate_names,
+            slice_size=slice_size,
+            progress=progress,
+        )
+        self._problem_fn = problem_fn
+        if self.io_api is None:
+            self.io_api = read_io_api_from_problem_file(problem_fn)
+        self._cache_model_sizes(model)
+
+    def solve(self, **run_kwargs: Any) -> Result:
+        """Run the prepared solver and return a :class:`Result`."""
+        if self.io_api == "direct" or self.solver_model is not None:
+            return self._run_direct(**run_kwargs)
+        if self._problem_fn is not None:
+            return self._run_file(**run_kwargs)
+        raise RuntimeError(
+            "Solver has not been built; call Solver.from_name(...) or _build() first."
+        )
+
+    def _run_direct(self, **run_kwargs: Any) -> Result:
+        """Run the pre-built native solver model. Override per-solver."""
+        raise NotImplementedError(
+            f"Direct API not implemented for {self.solver_name.value}"
+        )
+
+    def _run_file(self, **run_kwargs: Any) -> Result:
+        """Invoke the solver binary on ``self._problem_fn``. Override per-solver."""
+        raise NotImplementedError(
+            f"File-based API not implemented for {self.solver_name.value}"
+        )
 
     def solve_problem(
         self,
@@ -387,17 +563,19 @@ def solve_problem(
         env: EnvType | None = None,
         explicit_coordinate_names: bool = False,
     ) -> Result:
-        """
-        Solve a linear problem either from a model or a problem file.
-
-        Wraps around `self.solve_problem_from_model()` and
-        `self.solve_problem_from_file()` and calls the appropriate method
-        based on the input arguments (`model` or `problem_fn`).
-        """
+        """Deprecated. Use ``Solver.from_name(...).solve(...)`` or ``Model.solve(...)``."""
+        warnings.warn(
+            "Solver.solve_problem is deprecated and will be removed in a future "
+            "release. Use Solver.from_name(name, model, ...).solve(...) or "
+            "Model.solve(...) instead.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
         if problem_fn is not None and model is not None:
-            msg = "Both problem file and model are given. Please specify only one."
-            raise ValueError(msg)
-        elif model is not None:
+            raise ValueError(
+                "Both problem file and model are given. Please specify only one."
+            )
+        if model is not None:
             return self.solve_problem_from_model(
                 model=model,
                 solution_fn=solution_fn,
@@ -407,7 +585,7 @@ def solve_problem(
                 env=env,
                 explicit_coordinate_names=explicit_coordinate_names,
             )
-        elif problem_fn is not None:
+        if problem_fn is not None:
             return self.solve_problem_from_file(
                 problem_fn=problem_fn,
                 solution_fn=solution_fn,
@@ -416,30 +594,7 @@ def solve_problem(
                 basis_fn=basis_fn,
                 env=env,
             )
-        else:
-            msg = "No problem file or model specified."
-            raise ValueError(msg)
-
-    @property
-    def solver_name(self) -> SolverName:
-        return SolverName[self.__class__.__name__]
-
-
-class CBC(Solver[None]):
-    """
-    Solver subclass for the CBC solver.
-
-    Attributes
-    ----------
-    **solver_options
-        options for the given solver
-    """
-
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
+        raise ValueError("No problem file or model specified.")
 
     def solve_problem_from_model(
         self,
@@ -448,12 +603,38 @@ def solve_problem_from_model(
         log_fn: Path | None = None,
         warmstart_fn: Path | None = None,
         basis_fn: Path | None = None,
-        env: None = None,
+        env: EnvType | None = None,
         explicit_coordinate_names: bool = False,
         set_names: bool = True,
     ) -> Result:
-        msg = "Direct API not implemented for CBC"
-        raise NotImplementedError(msg)
+        """Deprecated shim that builds via ``_build_direct`` and runs via ``_run_direct``."""
+        warnings.warn(
+            "Solver.solve_problem_from_model is deprecated and will be removed in a "
+            "future release. Use Solver.from_name(name, model, io_api='direct', ...)"
+            ".solve(...) instead.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        if not self.supports(SolverFeature.DIRECT_API):
+            raise NotImplementedError(
+                f"Direct API not implemented for {self.solver_name.value}"
+            )
+        self.model = model
+        build_kwargs: dict[str, Any] = {
+            "explicit_coordinate_names": explicit_coordinate_names,
+            "set_names": set_names,
+            "log_fn": log_fn,
+        }
+        if env is not None:
+            build_kwargs["env"] = env
+        self._build_direct(**build_kwargs)
+        return self._run_direct(
+            solution_fn=solution_fn,
+            log_fn=log_fn,
+            warmstart_fn=warmstart_fn,
+            basis_fn=basis_fn,
+            env=env,
+        )
 
     def solve_problem_from_file(
         self,
@@ -462,34 +643,143 @@ def solve_problem_from_file(
         log_fn: Path | None = None,
         warmstart_fn: Path | None = None,
         basis_fn: Path | None = None,
-        env: None = None,
+        env: EnvType | None = None,
+    ) -> Result:
+        """Deprecated shim that caches ``problem_fn`` and runs via ``_run_file``."""
+        warnings.warn(
+            "Solver.solve_problem_from_file is deprecated and will be removed in a "
+            "future release. Use Solver.from_name(name, model, problem_fn=..., ...)"
+            ".solve(...) instead.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        problem_fn = (
+            Path(problem_fn) if not isinstance(problem_fn, Path) else problem_fn
+        )
+        self._problem_fn = problem_fn
+        self.io_api = read_io_api_from_problem_file(problem_fn)
+        return self._run_file(
+            solution_fn=solution_fn,
+            log_fn=log_fn,
+            warmstart_fn=warmstart_fn,
+            basis_fn=basis_fn,
+            env=env,
+        )
+
+    def _cache_model_labels(self, model: Model) -> None:
+        """Cache vlabels/clabels and total label counts for label-indexed solutions."""
+        self._vlabels = model.variables.label_index.vlabels
+        self._clabels = model.constraints.label_index.clabels
+        self._n_vars = model._xCounter
+        self._n_cons = model._cCounter
+
+    def _cache_model_sizes(self, model: Model) -> None:
+        """Cache total label counts only (file-based solvers parse names)."""
+        self._n_vars = model._xCounter
+        self._n_cons = model._cCounter
+
+    def update_solver_model(self, model: Model, **kwargs: Any) -> None:
+        raise NotImplementedError
+
+    def close(self) -> None:
+        if self._env_stack is not None:
+            self._env_stack.close()
+        self.env = None
+        self.solver_model = None
+        self._env_stack = None
+
+    def __del__(self) -> None:
+        with contextlib.suppress(Exception):
+            self.close()
+
+    def __repr__(self) -> str:
+        status = self.status.status.value if self.status is not None else "unsolved"
+        parts = [f"name={self.solver_name.value!r}", f"status={status!r}"]
+        if self.io_api is not None:
+            parts.append(f"io_api={self.io_api!r}")
+        if self.solver_model is not None:
+            parts.append("solver_model=loaded")
+        if self.env is not None:
+            parts.append("env=active")
+        if self.solution is not None:
+            parts.append(f"objective={self.solution.objective:.4g}")
+        if self.report is not None and self.report.runtime is not None:
+            parts.append(f"runtime={self.report.runtime:.3g}s")
+        return f"{type(self).__name__}({', '.join(parts)})"
+
+    def _make_result(
+        self,
+        status: Status,
+        solution: Solution | None,
+        solver_model: Any = None,
+        report: SolverReport | None = None,
     ) -> Result:
+        self.status = status
+        self.solution = solution
+        self.report = report
+        if solver_model is not None:
+            self.solver_model = solver_model
+        return Result(
+            status=status,
+            solution=solution,
+            solver_model=solver_model,
+            solver_name=self.solver_name.value,
+            report=report,
+        )
+
+    def safe_get_solution(
+        self, status: Status, func: Callable[[], Solution]
+    ) -> Solution:
+        """
+        Get solution from function call, if status is unknown still try to run it.
         """
-        Solve a linear problem from a problem file using the CBC solver.
+        if status.is_ok:
+            return func()
+        elif status.status == SolverStatus.unknown:
+            try:
+                logger.warning("Solution status unknown. Trying to parse solution.")
+                sol = func()
+                status.status = SolverStatus.ok
+                logger.warning("Solution parsed successfully.")
+                return sol
+            except Exception as e:
+                logger.error(f"Failed to parse solution: {e}")
+        return Solution()
 
-        The function reads the linear problem file and passes it to the solver.
-        If the solution is successful it returns variable solutions
-        and constraint dual values.
+    @property
+    def solver_name(self) -> SolverName:
+        return SolverName[self.__class__.__name__]
 
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path
-            Path to the solution file. This is necessary for solving with CBC.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
 
-        Returns
-        -------
-        Result
-        """
+class CBC(Solver[None]):
+    """
+    Solver subclass for the CBC solver.
+
+    Attributes
+    ----------
+    **solver_options
+        options for the given solver
+    """
+
+    display_name: ClassVar[str] = "CBC"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+        }
+    )
+
+    def _run_file(
+        self,
+        solution_fn: Path | None = None,
+        log_fn: Path | None = None,
+        warmstart_fn: Path | None = None,
+        basis_fn: Path | None = None,
+        env: None = None,
+        **kw: Any,
+    ) -> Result:
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         sense = read_sense_from_problem_file(problem_fn)
         io_api = read_io_api_from_problem_file(problem_fn)
 
@@ -588,9 +878,18 @@ def get_solver_solution() -> Solution:
             )
             variables_b = df.index.isin(variables)
 
-            sol = df[variables_b][2]
-            dual = df[~variables_b][3]
-
+            sol_df = df[variables_b]
+            dual_df = df[~variables_b]
+            sol = _solution_from_names(
+                sol_df[2].to_numpy(dtype=float),
+                sol_df.index.tolist(),
+                self._n_vars,
+            )
+            dual = _solution_from_names(
+                dual_df[3].to_numpy(dtype=float),
+                dual_df.index.tolist(),
+                self._n_cons,
+            )
             return Solution(sol, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
@@ -609,7 +908,13 @@ def get_solver_solution() -> Solution:
             runtime = float(m.group(1))
         CbcModel = namedtuple("CbcModel", ["mip_gap", "runtime"])
 
-        return Result(status, solution, CbcModel(mip_gap, runtime))
+        self.io_api = io_api
+        return self._make_result(
+            status,
+            solution,
+            solver_model=CbcModel(mip_gap, runtime),
+            report=SolverReport(runtime=runtime, mip_gap=mip_gap),
+        )
 
 
 class GLPK(Solver[None]):
@@ -622,65 +927,25 @@ class GLPK(Solver[None]):
         options for the given solver
     """
 
-    def __init(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
-
-    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: None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        msg = "Direct API not implemented for GLPK"
-        raise NotImplementedError(msg)
+    display_name: ClassVar[str] = "GLPK"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+        }
+    )
 
-    def solve_problem_from_file(
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the glpk solver.
-
-        This function reads the linear problem file and passes it to the
-        glpk solver. If the solution is successful it returns variable solutions
-        and constraint dual values.
-
-        For more information on the glpk solver options, see
-
-        https://kam.mff.cuni.cz/~elias/glpk.pdf
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path
-            Path to the solution file. This is necessary for solving with GLPK.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
-
-        Returns
-        -------
-        Result
-        """
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         CONDITION_MAP = {
             "integer optimal": "optimal",
             "integer undefined": "infeasible_or_unbounded",
@@ -740,7 +1005,8 @@ def solve_problem_from_file(
 
         if not os.path.exists(solution_fn):
             status = Status(SolverStatus.warning, TerminationCondition.unknown)
-            return Result(status, Solution())
+            self.io_api = io_api
+            return self._make_result(status, Solution())
 
         f = open(solution_fn)
 
@@ -764,23 +1030,31 @@ def get_solver_solution() -> Solution:
             dual_io = io.StringIO("".join(read_until_break(f))[:-2])
             dual_ = pd.read_fwf(dual_io)[1:].set_index("Row name")
             if "Marginal" in dual_:
-                dual = pd.to_numeric(dual_["Marginal"], "coerce").fillna(0)
+                dual = _solution_from_names(
+                    pd.to_numeric(dual_["Marginal"], "coerce")
+                    .fillna(0)
+                    .to_numpy(dtype=float),
+                    dual_.index.tolist(),
+                    self._n_cons,
+                )
             else:
                 logger.warning("Dual values of MILP couldn't be parsed")
-                dual = pd.Series(dtype=float)
+                dual = np.array([], dtype=float)
 
             sol_io = io.StringIO("".join(read_until_break(f))[:-2])
-            sol = (
-                pd.read_fwf(sol_io)[1:]
-                .set_index("Column name")["Activity"]
-                .astype(float)
+            sol_df = pd.read_fwf(sol_io)[1:].set_index("Column name")
+            sol = _solution_from_names(
+                sol_df["Activity"].astype(float).to_numpy(),
+                sol_df.index.tolist(),
+                self._n_vars,
             )
             f.close()
             return Solution(sol, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = maybe_adjust_objective_sign(solution, io_api, sense)
-        return Result(status, solution)
+        self.io_api = io_api
+        return self._make_result(status, solution)
 
 
 class Highs(Solver[None]):
@@ -803,54 +1077,29 @@ class Highs(Solver[None]):
         options for the given solver
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
+    display_name: ClassVar[str] = "HiGHS"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.DIRECT_API,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+            SolverFeature.SEMI_CONTINUOUS_VARIABLES,
+            SolverFeature.MIP_DUAL_BOUND_REPORT,
+        }
+    )
 
-    def solve_problem_from_model(
+    def _build_direct(
         self,
-        model: Model,
-        solution_fn: Path | None = None,
+        explicit_coordinate_names: bool = False,
+        set_names: bool = True,
         log_fn: Path | None = None,
-        warmstart_fn: Path | None = None,
-        basis_fn: Path | None = None,
-        env: None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        """
-        Solve a linear problem directly from a linopy model using the HiGHS solver.
-        Reads a linear problem file and passes it to the HiGHS solver.
-        If the solution is feasible the function returns the
-        objective, solution and dual constraint variables.
-
-        Parameters
-        ----------
-        model : linopy.model
-            Linopy model for the problem.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
-        explicit_coordinate_names : bool, optional
-            Transfer variable and constraint names to the solver (default: False)
-        set_names : bool, optional
-            Whether to set variable and constraint names (default: True).
-            Setting to False can significantly speed up model export.
-
-        Returns
-        -------
-        Result
-        """
-        # check for HiGHS solver compatibility
+        **kwargs: Any,
+    ) -> None:
+        model = self.model
+        assert model is not None
         if self.solver_options.get("solver") in [
             "simplex",
             "ipm",
@@ -864,70 +1113,130 @@ def solve_problem_from_model(
                 "Drop the solver option or use 'choose' to enable quadratic terms / integrality."
             )
 
-        h = model.to_highspy(
+        h = self._build_solver_model(
+            model,
             explicit_coordinate_names=explicit_coordinate_names,
             set_names=set_names,
         )
         self._set_solver_params(h, log_fn)
+        self.solver_model = h
+        self.io_api = "direct"
+        self.sense = model.sense
+        self._cache_model_labels(model)
+
+    @staticmethod
+    def _build_solver_model(
+        model: Model,
+        explicit_coordinate_names: bool = False,
+        set_names: bool = True,
+    ) -> highspy.Highs:
+        """Build a highspy.Highs instance that mirrors the linopy `model`."""
+        if model.variables.sos:
+            raise NotImplementedError(
+                "SOS constraints are not supported by the HiGHS direct API. "
+                "Use io_api='lp' instead."
+            )
+
+        M = model.matrices
+        h = highspy.Highs()
+        h.addVars(len(M.vlabels), M.lb, M.ub)
+        if (
+            len(model.binaries)
+            + len(model.integers)
+            + len(list(model.variables.semi_continuous))
+        ):
+            vtypes = M.vtypes
+            integrality_map = {"C": 0, "B": 1, "I": 1, "S": 2}
+            int_mask = (vtypes == "B") | (vtypes == "I") | (vtypes == "S")
+            labels = np.arange(len(vtypes))[int_mask]
+            integrality = np.array(
+                [integrality_map[v] for v in vtypes[int_mask]], dtype=np.int32
+            )
+            h.changeColsIntegrality(len(labels), labels, integrality)
+            if len(model.binaries):
+                labels = np.arange(len(vtypes))[vtypes == "B"]
+                n = len(labels)
+                h.changeColsBounds(
+                    n, labels, np.zeros_like(labels), np.ones_like(labels)
+                )
+
+        c = M.c
+        h.changeColsCost(len(c), np.arange(len(c), dtype=np.int32), c)
+
+        A = M.A
+        if A is not None:
+            A = A.tocsr()
+            num_cons = A.shape[0]
+            lower = np.where(M.sense != "<", M.b, -np.inf)
+            upper = np.where(M.sense != ">", M.b, np.inf)
+            h.addRows(num_cons, lower, upper, A.nnz, A.indptr, A.indices, A.data)
+
+        if set_names:
+            print_variables, print_constraints = linopy.io.get_printers_scalar(
+                model, explicit_coordinate_names=explicit_coordinate_names
+            )
+            lp = h.getLp()
+            lp.col_names_ = print_variables(M.vlabels)
+            if len(M.clabels):
+                lp.row_names_ = print_constraints(M.clabels)
+            h.passModel(lp)
+
+        Q = M.Q
+        if Q is not None:
+            Q = triu(Q).tocsr()
+            num_vars = Q.shape[0]
+            h.passHessian(num_vars, Q.nnz, 1, Q.indptr, Q.indices, Q.data)
+
+        if model.objective.sense == "max":
+            h.changeObjectiveSense(highspy.ObjSense.kMaximize)
+
+        return h
 
+    def _run_direct(
+        self,
+        solution_fn: Path | None = None,
+        log_fn: Path | None = None,
+        warmstart_fn: Path | None = None,
+        basis_fn: Path | None = None,
+        env: Any = None,
+        **kw: Any,
+    ) -> Result:
         return self._solve(
-            h,
-            solution_fn,
-            warmstart_fn,
-            basis_fn,
-            model=model,
-            io_api="direct",
-            sense=model.sense,
+            self.solver_model,
+            solution_fn=solution_fn,
+            warmstart_fn=warmstart_fn,
+            basis_fn=basis_fn,
+            io_api=self.io_api,
+            sense=self.sense,
         )
 
-    def solve_problem_from_file(
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the HiGHS solver.
-        Reads a linear problem file and passes it to the HiGHS solver.
-        If the solution is feasible the function returns the
-        objective, solution and dual constraint variables.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
-
-        Returns
-        -------
-        Result
-        """
-
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         problem_fn_ = path_to_string(problem_fn)
         h = highspy.Highs()
         self._set_solver_params(h, log_fn)
 
         h.readModel(problem_fn_)
+        self.solver_model = h
+        self.io_api = read_io_api_from_problem_file(problem_fn)
 
         return self._solve(
             h,
             solution_fn,
             warmstart_fn,
             basis_fn,
-            io_api=read_io_api_from_problem_file(problem_fn),
+            io_api=self.io_api,
             sense=read_sense_from_problem_file(problem_fn),
+            from_file=True,
         )
 
     def _set_solver_params(
@@ -948,9 +1257,9 @@ def _solve(
         solution_fn: Path | None = None,
         warmstart_fn: Path | None = None,
         basis_fn: Path | None = None,
-        model: Model | None = None,
         io_api: str | None = None,
         sense: str | None = None,
+        from_file: bool = False,
     ) -> Result:
         """
         Solve a linear problem from a HiGHS object.
@@ -968,12 +1277,13 @@ def _solve(
             Path to the warmstart file.
         basis_fn : Path, optional
             Path to the basis file.
-        model : linopy.model, optional
-            Linopy model for the problem.
         io_api: str
             io_api of the problem. For direct API from linopy model this is "direct".
         sense: str
             "min" or "max"
+        from_file: bool
+            ``True`` when ``h`` was populated via ``readModel`` — HiGHS may have
+            reordered columns/rows, so values are re-permuted using parsed names.
 
         Returns
         -------
@@ -1025,20 +1335,39 @@ def _solve(
         def get_solver_solution() -> Solution:
             objective = h.getObjectiveValue()
             solution = h.getSolution()
-
-            if model is not None:
-                sol = pd.Series(solution.col_value, model.matrices.vlabels, dtype=float)
-                dual = pd.Series(solution.row_dual, model.matrices.clabels, dtype=float)
+            sol = np.asarray(solution.col_value, dtype=float)
+            dual = np.asarray(solution.row_dual, dtype=float)
+            if from_file:
+                lp = h.getLp()
+                sol = _solution_from_names(sol, lp.col_names_, self._n_vars)
+                dual = _solution_from_names(dual, lp.row_names_, self._n_cons)
             else:
-                sol = pd.Series(solution.col_value, h.getLp().col_names_, dtype=float)
-                dual = pd.Series(solution.row_dual, h.getLp().row_names_, dtype=float)
-
+                sol = _solution_from_labels(sol, self._vlabels, self._n_vars)
+                dual = _solution_from_labels(dual, self._clabels, self._n_cons)
             return Solution(sol, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = maybe_adjust_objective_sign(solution, io_api, sense)
 
-        return Result(status, solution, h)
+        runtime: float | None = None
+        mip_gap: float | None = None
+        dual_bound: float | None = None
+        with contextlib.suppress(Exception):
+            runtime = float(h.getRunTime())
+        with contextlib.suppress(Exception):
+            mip_gap = float(h.getInfo().mip_gap)
+        with contextlib.suppress(Exception):
+            dual_bound = float(h.getInfo().mip_dual_bound)
+
+        self.io_api = io_api
+        return self._make_result(
+            status,
+            solution,
+            solver_model=h,
+            report=SolverReport(
+                runtime=runtime, mip_gap=mip_gap, dual_bound=dual_bound
+            ),
+        )
 
 
 class Gurobi(Solver["gurobipy.Env | dict[str, Any] | None"]):
@@ -1051,132 +1380,169 @@ class Gurobi(Solver["gurobipy.Env | dict[str, Any] | None"]):
         options for the given solver
     """
 
-    def __init__(
+    display_name: ClassVar[str] = "Gurobi"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.DIRECT_API,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+            SolverFeature.IIS_COMPUTATION,
+            SolverFeature.SOS_CONSTRAINTS,
+            SolverFeature.SEMI_CONTINUOUS_VARIABLES,
+            SolverFeature.SOLVER_ATTRIBUTE_ACCESS,
+            SolverFeature.MIP_DUAL_BOUND_REPORT,
+        }
+    )
+
+    def _resolve_env(self, env: gurobipy.Env | dict[str, Any] | None) -> gurobipy.Env:
+        self.close()
+        self._env_stack = contextlib.ExitStack()
+        if env is None:
+            resolved = self._env_stack.enter_context(gurobipy.Env())
+        elif isinstance(env, dict):
+            resolved = self._env_stack.enter_context(gurobipy.Env(params=env))
+        else:
+            resolved = env
+        self.env = resolved
+        return resolved
+
+    def _build_direct(
         self,
-        **solver_options: Any,
+        explicit_coordinate_names: bool = False,
+        env: gurobipy.Env | dict[str, Any] | None = None,
+        set_names: bool = True,
+        **kwargs: Any,
     ) -> None:
-        super().__init__(**solver_options)
+        model = self.model
+        assert model is not None
+        env_ = self._resolve_env(env)
+        m = self._build_solver_model(
+            model,
+            env=env_,
+            explicit_coordinate_names=explicit_coordinate_names,
+            set_names=set_names,
+        )
+        self.solver_model = m
+        self.io_api = "direct"
+        self.sense = model.sense
+        self._cache_model_labels(model)
 
-    def solve_problem_from_model(
-        self,
+    @staticmethod
+    def _build_solver_model(
         model: Model,
-        solution_fn: Path | None = None,
-        log_fn: Path | None = None,
-        warmstart_fn: Path | None = None,
-        basis_fn: Path | None = None,
-        env: gurobipy.Env | dict[str, Any] | None = None,
+        env: gurobipy.Env | None = None,
         explicit_coordinate_names: bool = False,
         set_names: bool = True,
-    ) -> Result:
-        """
-        Solve a linear problem directly from a linopy model using the Gurobi solver.
-        Reads a problem file and passes it to the Gurobi solver.
-        This function communicates with gurobi using the gurobipy package.
+    ) -> gurobipy.Model:
+        """Build a gurobipy.Model that mirrors the linopy `model`."""
+        model.constraints.sanitize_missings()
+        gm = gurobipy.Model(env=env)
 
-        Parameters
-        ----------
-        model : linopy.model
-            Linopy model for the problem.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : gurobipy.Env or dict, optional
-            Gurobi environment for the solver, pass env directly or kwargs for creation.
-        explicit_coordinate_names : bool, optional
-            Transfer variable and constraint names to the solver (default: False)
-        set_names : bool, optional
-            Whether to set variable and constraint names (default: True).
-            Setting to False can significantly speed up model export.
+        M = model.matrices
 
-        Returns
-        -------
-        Result
-        """
-        with contextlib.ExitStack() as stack:
-            if env is None:
-                env_ = stack.enter_context(gurobipy.Env())
-            elif isinstance(env, dict):
-                env_ = stack.enter_context(gurobipy.Env(params=env))
-            else:
-                env_ = env
-
-            m = model.to_gurobipy(
-                env=env_,
-                explicit_coordinate_names=explicit_coordinate_names,
-                set_names=set_names,
+        kwargs: dict[str, Any] = {}
+        if set_names:
+            print_variables, print_constraints = linopy.io.get_printers_scalar(
+                model, explicit_coordinate_names=explicit_coordinate_names
             )
+            kwargs["name"] = print_variables(M.vlabels)
+        if (
+            len(model.binaries.labels)
+            + len(model.integers.labels)
+            + len(list(model.variables.semi_continuous))
+        ):
+            kwargs["vtype"] = M.vtypes
+        x = gm.addMVar(M.vlabels.shape, M.lb, M.ub, **kwargs)
+
+        if model.is_quadratic:
+            gm.setObjective(0.5 * x.T @ M.Q @ x + M.c @ x)
+        else:
+            gm.setObjective(M.c @ x)
+
+        if model.objective.sense == "max":
+            gm.ModelSense = -1
+
+        if len(model.constraints):
+            c = gm.addMConstr(M.A, x, M.sense, M.b)
+            if set_names:
+                names = print_constraints(M.clabels)
+                c.setAttr("ConstrName", names)
+
+        if model.variables.sos:
+            for var_name in model.variables.sos:
+                var = model.variables.sos[var_name]
+                sos_type: int = var.attrs[SOS_TYPE_ATTR]
+                sos_dim: str = var.attrs[SOS_DIM_ATTR]
+
+                def add_sos(s: xr.DataArray, sos_type: int, sos_dim: str) -> None:
+                    s = s.squeeze()
+                    indices = s.values.flatten().tolist()
+                    weights = s.coords[sos_dim].values.tolist()
+                    gm.addSOS(sos_type, x[indices].tolist(), weights)
+
+                others = [dim for dim in var.labels.dims if dim != sos_dim]
+                if not others:
+                    add_sos(var.labels, sos_type, sos_dim)
+                else:
+                    stacked = var.labels.stack(_sos_group=others)
+                    for _, s in stacked.groupby("_sos_group"):
+                        add_sos(s.unstack("_sos_group"), sos_type, sos_dim)
 
-            return self._solve(
-                m,
-                solution_fn=solution_fn,
-                log_fn=log_fn,
-                warmstart_fn=warmstart_fn,
-                basis_fn=basis_fn,
-                io_api="direct",
-                sense=model.sense,
-            )
+        gm.update()
+        return gm
 
-    def solve_problem_from_file(
+    def _run_direct(
         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: gurobipy.Env | dict[str, Any] | None = None,
+        env: Any = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the Gurobi solver.
-        Reads a problem file and passes it to the Gurobi solver.
-        This function communicates with gurobi using the gurobipy package.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : gurobipy.Env or dict, optional
-            Gurobi environment for the solver, pass env directly or kwargs for creation.
+        return self._solve(
+            self.solver_model,
+            solution_fn=solution_fn,
+            log_fn=log_fn,
+            warmstart_fn=warmstart_fn,
+            basis_fn=basis_fn,
+            io_api=self.io_api,
+            sense=self.sense,
+        )
 
-        Returns
-        -------
-        Result
-        """
+    def _run_file(
+        self,
+        solution_fn: Path | None = None,
+        log_fn: Path | None = None,
+        warmstart_fn: Path | None = None,
+        basis_fn: Path | None = None,
+        env: gurobipy.Env | dict[str, Any] | None = None,
+        **kw: Any,
+    ) -> Result:
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         sense = read_sense_from_problem_file(problem_fn)
         io_api = read_io_api_from_problem_file(problem_fn)
         problem_fn_ = path_to_string(problem_fn)
 
-        with contextlib.ExitStack() as stack:
-            if env is None:
-                env_ = stack.enter_context(gurobipy.Env())
-            elif isinstance(env, dict):
-                env_ = stack.enter_context(gurobipy.Env(params=env))
-            else:
-                env_ = env
-
-            m = gurobipy.read(problem_fn_, env=env_)
+        env_ = self._resolve_env(env)
+        m = gurobipy.read(problem_fn_, env=env_)
+        self.solver_model = m
+        self.io_api = io_api
 
-            return self._solve(
-                m,
-                solution_fn=solution_fn,
-                log_fn=log_fn,
-                warmstart_fn=warmstart_fn,
-                basis_fn=basis_fn,
-                io_api=io_api,
-                sense=sense,
-            )
+        return self._solve(
+            m,
+            solution_fn=solution_fn,
+            log_fn=log_fn,
+            warmstart_fn=warmstart_fn,
+            basis_fn=basis_fn,
+            io_api=io_api,
+            sense=sense,
+            from_file=True,
+        )
 
     def _solve(
         self,
@@ -1187,6 +1553,7 @@ def _solve(
         basis_fn: Path | None,
         io_api: str | None,
         sense: str | None,
+        from_file: bool = False,
     ) -> Result:
         """
         Solve a linear problem from a Gurobi object.
@@ -1264,22 +1631,54 @@ def _solve(
         def get_solver_solution() -> Solution:
             objective = m.ObjVal
 
-            sol = pd.Series({v.VarName: v.x for v in m.getVars()}, dtype=float)  # type: ignore
+            vars_ = m.getVars()
+            sol = np.array([v.X for v in vars_], dtype=float)
+            if from_file:
+                sol = _solution_from_names(
+                    sol, [v.VarName for v in vars_], self._n_vars
+                )
+            else:
+                sol = _solution_from_labels(sol, self._vlabels, self._n_vars)
 
             try:
-                dual = pd.Series(
-                    {c.ConstrName: c.Pi for c in m.getConstrs()}, dtype=float
-                )
+                constrs = m.getConstrs()
+                dual = np.array([c.Pi for c in constrs], dtype=float)
+                if from_file:
+                    dual = _solution_from_names(
+                        dual,
+                        [c.ConstrName for c in constrs],
+                        self._n_cons,
+                    )
+                else:
+                    dual = _solution_from_labels(dual, self._clabels, self._n_cons)
             except AttributeError:
                 logger.warning("Dual values of MILP couldn't be parsed")
-                dual = pd.Series(dtype=float)
+                dual = np.array([], dtype=float)
 
             return Solution(sol, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = maybe_adjust_objective_sign(solution, io_api, sense)
 
-        return Result(status, solution, m)
+        runtime: float | None = None
+        mip_gap: float | None = None
+        dual_bound: float | None = None
+        with contextlib.suppress(Exception):
+            runtime = float(m.Runtime)
+        with contextlib.suppress(Exception):
+            mip_gap = float(m.MIPGap)
+        with contextlib.suppress(Exception):
+            dual_bound = float(m.ObjBound)
+
+        self.io_api = io_api
+        return self._make_result(
+            status,
+            solution,
+            solver_model=m,
+            report=SolverReport(
+                runtime=runtime, mip_gap=mip_gap, dual_bound=dual_bound
+            ),
+        )
 
 
 class Cplex(Solver[None]):
@@ -1296,61 +1695,29 @@ class Cplex(Solver[None]):
         options for the given solver
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
-
-    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: None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        msg = "Direct API not implemented for Cplex"
-        raise NotImplementedError(msg)
+    display_name: ClassVar[str] = "CPLEX"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOS_CONSTRAINTS,
+            SolverFeature.SEMI_CONTINUOUS_VARIABLES,
+        }
+    )
 
-    def solve_problem_from_file(
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the cplex solver.
-
-        This function reads the linear problem file and passes it to the cplex
-        solver. If the solution is successful it returns variable solutions and
-        constraint dual values. Cplex must be installed for using this function.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
-
-        Returns
-        -------
-        Result
-        """
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         CONDITION_MAP = {
             "integer optimal solution": "optimal",
             "integer optimal, tolerance": "optimal",
@@ -1418,27 +1785,30 @@ def get_solver_solution() -> Solution:
 
             objective = m.solution.get_objective_value()
 
-            solution = pd.Series(
-                m.solution.get_values(), m.variables.get_names(), dtype=float
+            solution = _solution_from_names(
+                np.asarray(m.solution.get_values(), dtype=float),
+                m.variables.get_names(),
+                self._n_vars,
             )
 
             try:
-                dual = pd.Series(
-                    m.solution.get_dual_values(),
+                dual = _solution_from_names(
+                    np.asarray(m.solution.get_dual_values(), dtype=float),
                     m.linear_constraints.get_names(),
-                    dtype=float,
+                    self._n_cons,
                 )
             except Exception:
                 logger.warning(
                     "Dual values not available (e.g. barrier solution without crossover)"
                 )
-                dual = pd.Series(dtype=float)
+                dual = np.array([], dtype=float)
             return Solution(solution, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = maybe_adjust_objective_sign(solution, io_api, sense)
 
-        return Result(status, solution, m)
+        self.io_api = io_api
+        return self._make_result(status, solution, solver_model=m)
 
 
 class SCIP(Solver[None]):
@@ -1451,59 +1821,28 @@ class SCIP(Solver[None]):
         options for the given solver
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
-
-    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: None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        msg = "Direct API not implemented for SCIP"
-        raise NotImplementedError(msg)
+    display_name: ClassVar[str] = "SCIP"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+        }
+    )
 
-    def solve_problem_from_file(
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the scip solver.
-
-        This function communicates with scip using the pyscipopt package.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
-
-        Returns
-        -------
-        Result
-        """
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         CONDITION_MAP: dict[str, TerminationCondition] = {
             # https://github.com/scipopt/scip/blob/b2bac412222296ff2b7f2347bb77d5fc4e05a2a1/src/scip/type_stat.h#L40
             "inforunbd": TerminationCondition.infeasible_or_unbounded,
@@ -1570,29 +1909,32 @@ def get_solver_solution() -> Solution:
             vars_to_ignore = {"quadobjvar", "qmatrixvar", "quadobj", "qmatrix"}
 
             s = m.getSols()[0]
-            sol = pd.Series(
-                {v.name: s[v] for v in m.getVars() if v.name not in vars_to_ignore}
+            kept_vars = [v for v in m.getVars() if v.name not in vars_to_ignore]
+            sol = _solution_from_names(
+                np.array([s[v] for v in kept_vars], dtype=float),
+                [v.name for v in kept_vars],
+                self._n_vars,
             )
 
             cons = m.getConss(False)
             if len(cons) != 0:
-                dual = pd.Series(
-                    {
-                        c.name: m.getDualSolVal(c)
-                        for c in cons
-                        if c.name not in vars_to_ignore
-                    }
+                kept_cons = [c for c in cons if c.name not in vars_to_ignore]
+                dual = _solution_from_names(
+                    np.array([m.getDualSolVal(c) for c in kept_cons], dtype=float),
+                    [c.name for c in kept_cons],
+                    self._n_cons,
                 )
             else:
                 logger.warning("Dual values not available (is this an MILP?)")
-                dual = pd.Series(dtype=float)
+                dual = np.array([], dtype=float)
 
             return Solution(sol, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = maybe_adjust_objective_sign(solution, io_api, sense)
 
-        return Result(status, solution, m)
+        self.io_api = io_api
+        return self._make_result(status, solution, solver_model=m)
 
 
 class Xpress(Solver[None]):
@@ -1608,62 +1950,35 @@ class Xpress(Solver[None]):
         options for the given solver
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
+    display_name: ClassVar[str] = "FICO Xpress"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+            SolverFeature.IIS_COMPUTATION,
+        }
+    )
 
-    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: None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        msg = "Direct API not implemented for Xpress"
-        raise NotImplementedError(msg)
+    @classmethod
+    def runtime_features(cls) -> frozenset[SolverFeature]:
+        if _installed_version_in("xpress", ">=9.8.0"):
+            return frozenset({SolverFeature.GPU_ACCELERATION})
+        return frozenset()
 
-    def solve_problem_from_file(
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the Xpress solver.
-
-        This function reads the linear problem file and passes it to
-        the Xpress solver. If the solution is successful it returns
-        variable solutions and constraint dual values. The `xpress` module
-        must be installed for using this function.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
-
-        Returns
-        -------
-        Result
-        """
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         CONDITION_MAP = {
             xpress.SolStatus.NOTFOUND: "unknown",
             xpress.SolStatus.OPTIMAL: "optimal",
@@ -1733,40 +2048,36 @@ def solve_problem_from_file(
         def get_solver_solution() -> Solution:
             objective = m.attributes.objval
 
-            try:  # Try new API first
-                var = m.getNameList(xpress_Namespaces.COLUMN, 0, m.attributes.cols - 1)
-            except AttributeError:  # Fallback to old API
-                var = m.getnamelist(xpress_Namespaces.COLUMN, 0, m.attributes.cols - 1)
-            sol = pd.Series(m.getSolution(), index=var, dtype=float)
+            sol = _solution_from_names(
+                np.asarray(m.getSolution(), dtype=float),
+                [v.name for v in m.getVariable()],
+                self._n_vars,
+            )
 
             try:
                 if m.attributes.rows == 0:
-                    dual = pd.Series(dtype=float)
+                    dual = np.array([], dtype=float)
                 else:
                     try:  # Try new API first
                         _dual = m.getDuals()
                     except AttributeError:  # Fallback to old API
                         _dual = m.getDual()
-
-                    try:  # Try new API first
-                        constraints = m.getNameList(
-                            xpress_Namespaces.ROW, 0, m.attributes.rows - 1
-                        )
-                    except AttributeError:  # Fallback to old API
-                        constraints = m.getnamelist(
-                            xpress_Namespaces.ROW, 0, m.attributes.rows - 1
-                        )
-                    dual = pd.Series(_dual, index=constraints, dtype=float)
+                    dual = _solution_from_names(
+                        np.asarray(_dual, dtype=float),
+                        [c.name for c in m.getConstraint()],
+                        self._n_cons,
+                    )
             except (xpress.SolverError, xpress.ModelError, SystemError):
                 logger.warning("Dual values of MILP couldn't be parsed")
-                dual = pd.Series(dtype=float)
+                dual = np.array([], dtype=float)
 
             return Solution(sol, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = maybe_adjust_objective_sign(solution, io_api, sense)
 
-        return Result(status, solution, m)
+        self.io_api = io_api
+        return self._make_result(status, solution, solver_model=m)
 
 
 KnitroResult = namedtuple(
@@ -1788,25 +2099,17 @@ class Knitro(Solver[None]):
         options for the given solver
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
-
-    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: None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        msg = "Direct API not implemented for Knitro"
-        raise NotImplementedError(msg)
+    display_name: ClassVar[str] = "Artelys Knitro"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+            SolverFeature.MIP_DUAL_BOUND_REPORT,
+        }
+    )
 
     @staticmethod
     def _set_option(kc: Any, name: str, value: Any) -> None:
@@ -1830,11 +2133,10 @@ def _extract_values(
         kc: Any,
         get_count_fn: Callable[..., Any],
         get_values_fn: Callable[..., Any],
-        get_names_fn: Callable[..., Any],
-    ) -> pd.Series:
+    ) -> np.ndarray:
         n = int(get_count_fn(kc))
         if n == 0:
-            return pd.Series(dtype=float)
+            return np.array([], dtype=float)
 
         try:
             # Compatible with KNITRO >= 15
@@ -1843,40 +2145,19 @@ def _extract_values(
             # Fallback for older wrappers requiring explicit indices
             values = get_values_fn(kc, list(range(n)))
 
-        names = list(get_names_fn(kc))
-        return pd.Series(values, index=names, dtype=float)
+        return np.asarray(values, dtype=float)
 
-    def solve_problem_from_file(
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the Knitro solver.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver.
-
-        Returns
-        -------
-        Result
-        """
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         CONDITION_MAP: dict[int, TerminationCondition] = {
             0: TerminationCondition.optimal,
             -100: TerminationCondition.suboptimal,
@@ -1971,19 +2252,23 @@ def get_solver_solution() -> Solution:
                     kc,
                     knitro.KN_get_number_vars,
                     knitro.KN_get_var_primal_values,
-                    knitro.KN_get_var_names,
                 )
+                n_vars = int(knitro.KN_get_number_vars(kc))
+                var_names = [knitro.KN_get_var_names(kc, i) for i in range(n_vars)]
+                sol = _solution_from_names(sol, var_names, self._n_vars)
 
                 try:
                     dual = self._extract_values(
                         kc,
                         knitro.KN_get_number_cons,
                         knitro.KN_get_con_dual_values,
-                        knitro.KN_get_con_names,
                     )
+                    n_cons = int(knitro.KN_get_number_cons(kc))
+                    con_names = [knitro.KN_get_con_names(kc, i) for i in range(n_cons)]
+                    dual = _solution_from_names(dual, con_names, self._n_cons)
                 except Exception:
                     logger.warning("Dual values couldn't be parsed")
-                    dual = pd.Series(dtype=float)
+                    dual = np.array([], dtype=float)
 
                 return Solution(sol, dual, objective)
 
@@ -1994,25 +2279,28 @@ def get_solver_solution() -> Solution:
                 solution_fn.parent.mkdir(exist_ok=True)
                 knitro.KN_write_mps_file(kc, path_to_string(solution_fn))
 
-            return Result(
+            knitro_model = KnitroResult(
+                reported_runtime=reported_runtime,
+                mip_relaxation_bnd=mip_relaxation_bnd,
+                mip_number_nodes=mip_number_nodes,
+                mip_number_solves=mip_number_solves,
+                mip_rel_gap=mip_rel_gap,
+                mip_abs_gap=mip_abs_gap,
+                abs_feas_error=abs_feas_error,
+                rel_feas_error=rel_feas_error,
+                abs_opt_error=abs_opt_error,
+                rel_opt_error=rel_opt_error,
+                n_vars=n_vars,
+                n_cons=n_cons,
+                n_integer_vars=n_integer_vars,
+                n_continuous_vars=n_continuous_vars,
+            )
+            self.io_api = io_api
+            return self._make_result(
                 status,
                 solution,
-                KnitroResult(
-                    reported_runtime=reported_runtime,
-                    mip_relaxation_bnd=mip_relaxation_bnd,
-                    mip_number_nodes=mip_number_nodes,
-                    mip_number_solves=mip_number_solves,
-                    mip_rel_gap=mip_rel_gap,
-                    mip_abs_gap=mip_abs_gap,
-                    abs_feas_error=abs_feas_error,
-                    rel_feas_error=rel_feas_error,
-                    abs_opt_error=abs_opt_error,
-                    rel_opt_error=rel_opt_error,
-                    n_vars=n_vars,
-                    n_cons=n_cons,
-                    n_integer_vars=n_integer_vars,
-                    n_continuous_vars=n_continuous_vars,
-                ),
+                solver_model=knitro_model,
+                report=SolverReport(runtime=reported_runtime, mip_gap=mip_rel_gap),
             )
         finally:
             with contextlib.suppress(Exception):
@@ -2042,135 +2330,180 @@ class Mosek(Solver[None]):
         options for the given solver
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
+    display_name: ClassVar[str] = "MOSEK"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.DIRECT_API,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+        }
+    )
 
-    def solve_problem_from_model(
+    def _run_direct(
         self,
-        model: Model,
         solution_fn: Path | None = None,
         log_fn: Path | None = None,
         warmstart_fn: Path | None = None,
         basis_fn: Path | None = None,
-        env: None = None,
+        env: Any = None,
+        **kw: Any,
+    ) -> Result:
+        return self._solve(
+            self.solver_model,
+            solution_fn=solution_fn,
+            log_fn=log_fn,
+            warmstart_fn=warmstart_fn,
+            basis_fn=basis_fn,
+            io_api=self.io_api,
+            sense=self.sense,
+        )
+
+    def _build_direct(
+        self,
         explicit_coordinate_names: bool = False,
         set_names: bool = True,
-    ) -> Result:
-        """
-        Solve a linear problem directly from a linopy model using the MOSEK solver.
+        **kwargs: Any,
+    ) -> None:
+        model = self.model
+        assert model is not None
+        self.close()
+        self._env_stack = contextlib.ExitStack()
+        task = self._env_stack.enter_context(mosek.Task())
+        m = self._build_solver_model(
+            model,
+            task,
+            explicit_coordinate_names=explicit_coordinate_names,
+            set_names=set_names,
+        )
+        self.solver_model = m
+        self.io_api = "direct"
+        self.sense = model.sense
+        self._cache_model_labels(model)
 
-        Parameters
-        ----------
-        model : linopy.model
-            Linopy model for the problem.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional, deprecated
-            Deprecated. This parameter is ignored. MOSEK now uses the global
-            environment automatically. Will be removed in a future version.
-        explicit_coordinate_names : bool, optional
-            Transfer variable and constraint names to the solver (default: False)
-        set_names : bool, optional
-            Whether to set variable and constraint names (default: True).
-            Setting to False can significantly speed up model export.
+    @staticmethod
+    def _build_solver_model(
+        model: Model,
+        task: mosek.Task,
+        explicit_coordinate_names: bool = False,
+        set_names: bool = True,
+    ) -> mosek.Task:
+        """Populate an empty MOSEK task with the contents of `model`."""
+        if model.variables.sos:
+            raise NotImplementedError("SOS constraints are not supported by MOSEK.")
+        if model.variables.semi_continuous:
+            raise NotImplementedError(
+                "Semi-continuous variables are not supported by MOSEK. "
+                "Use a solver that supports them (gurobi, cplex, highs)."
+            )
 
-        Returns
-        -------
-        Result
-        """
+        task.appendvars(model.nvars)
+        task.appendcons(model.ncons)
 
-        if env is not None:
-            warnings.warn(
-                "The 'env' parameter in solve_problem_from_model is deprecated and will be "
-                "removed in a future version. MOSEK now uses the global environment "
-                "automatically, avoiding unnecessary license checkouts.",
-                DeprecationWarning,
-                stacklevel=2,
+        M = model.matrices
+
+        if set_names:
+            print_variables, print_constraints = linopy.io.get_printers_scalar(
+                model, explicit_coordinate_names=explicit_coordinate_names
             )
-        with mosek.Task() as m:
-            m = model.to_mosek(
-                m,
-                explicit_coordinate_names=explicit_coordinate_names,
-                set_names=set_names,
+            labels = print_variables(M.vlabels)
+            task.generatevarnames(
+                np.arange(0, len(labels)), "%0", [len(labels)], None, [0], labels
             )
 
-            return self._solve(
-                m,
-                solution_fn=solution_fn,
-                log_fn=log_fn,
-                warmstart_fn=warmstart_fn,
-                basis_fn=basis_fn,
-                io_api="direct",
-                sense=model.sense,
+        bkx = [
+            (
+                (
+                    (mosek.boundkey.ra if lb < ub else mosek.boundkey.fx)
+                    if ub < np.inf
+                    else mosek.boundkey.lo
+                )
+                if (lb > -np.inf)
+                else (mosek.boundkey.up if (ub < np.inf) else mosek.boundkey.fr)
             )
+            for (lb, ub) in zip(M.lb, M.ub)
+        ]
+        blx = [b if b > -np.inf else 0.0 for b in M.lb]
+        bux = [b if b < np.inf else 0.0 for b in M.ub]
+        task.putvarboundslice(0, model.nvars, bkx, blx, bux)
+
+        if len(model.binaries.labels) + len(model.integers.labels) > 0:
+            idx = [i for (i, v) in enumerate(M.vtypes) if v in ["B", "I"]]
+            task.putvartypelist(idx, [mosek.variabletype.type_int] * len(idx))
+            if len(model.binaries.labels) > 0:
+                bidx = [i for (i, v) in enumerate(M.vtypes) if v == "B"]
+                task.putvarboundlistconst(bidx, mosek.boundkey.ra, 0.0, 1.0)
+
+        if len(model.constraints) > 0:
+            if set_names:
+                names = print_constraints(M.clabels)
+                for i, n in enumerate(names):
+                    task.putconname(i, n)
+            bkc = [
+                (
+                    (mosek.boundkey.up if b < np.inf else mosek.boundkey.fr)
+                    if s == "<"
+                    else (
+                        (mosek.boundkey.lo if b > -np.inf else mosek.boundkey.up)
+                        if s == ">"
+                        else mosek.boundkey.fx
+                    )
+                )
+                for s, b in zip(M.sense, M.b)
+            ]
+            blc = [b if b > -np.inf else 0.0 for b in M.b]
+            buc = [b if b < np.inf else 0.0 for b in M.b]
+            if M.A is not None:
+                A = M.A.tocsr()
+                task.putarowslice(
+                    0, model.ncons, A.indptr[:-1], A.indptr[1:], A.indices, A.data
+                )
+                task.putconboundslice(0, model.ncons, bkc, blc, buc)
 
-    def solve_problem_from_file(
+        if M.Q is not None:
+            Q = (0.5 * tril(M.Q + M.Q.transpose())).tocoo()
+            task.putqobj(Q.row, Q.col, Q.data)
+        task.putclist(list(np.arange(model.nvars)), M.c)
+
+        if model.objective.sense == "max":
+            task.putobjsense(mosek.objsense.maximize)
+        else:
+            task.putobjsense(mosek.objsense.minimize)
+        return task
+
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the MOSEK solver. Both mps and
-        lp files are supported; MPS does not support quadratic terms.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional, deprecated
-            Deprecated. This parameter is ignored. MOSEK now uses the global
-            environment automatically. Will be removed in a future version.
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
+        self.close()
+        self._env_stack = contextlib.ExitStack()
+        m = self._env_stack.enter_context(mosek.Task())
+        sense = read_sense_from_problem_file(problem_fn)
+        io_api = read_io_api_from_problem_file(problem_fn)
+        problem_fn_ = path_to_string(problem_fn)
+        m.readdata(problem_fn_)
+        self.solver_model = m
+        self.io_api = io_api
 
-        Returns
-        -------
-        Result
-        """
-        if env is not None:
-            warnings.warn(
-                "The 'env' parameter in solve_problem_from_file is deprecated and will be "
-                "removed in a future version. MOSEK now uses the global environment "
-                "automatically, avoiding unnecessary license checkouts.",
-                DeprecationWarning,
-                stacklevel=2,
-            )
-        with mosek.Task() as m:
-            # read sense and io_api from problem file
-            sense = read_sense_from_problem_file(problem_fn)
-            io_api = read_io_api_from_problem_file(problem_fn)
-            # for Mosek solver, the path needs to be a string
-            problem_fn_ = path_to_string(problem_fn)
-            m.readdata(problem_fn_)
-
-            return self._solve(
-                m,
-                solution_fn=solution_fn,
-                log_fn=log_fn,
-                warmstart_fn=warmstart_fn,
-                basis_fn=basis_fn,
-                io_api=io_api,
-                sense=sense,
-            )
+        return self._solve(
+            m,
+            solution_fn=solution_fn,
+            log_fn=log_fn,
+            warmstart_fn=warmstart_fn,
+            basis_fn=basis_fn,
+            io_api=io_api,
+            sense=sense,
+            from_file=True,
+        )
 
     def _solve(
         self,
@@ -2181,6 +2514,7 @@ def _solve(
         basis_fn: Path | None,
         io_api: str | None,
         sense: str | None,
+        from_file: bool = False,
     ) -> Result:
         """
         Solve a linear problem from a Mosek task object.
@@ -2349,24 +2683,41 @@ def _solve(
         def get_solver_solution() -> Solution:
             objective = m.getprimalobj(soltype)
 
-            sol = m.getxx(soltype)
-            sol = {m.getvarname(i): sol[i] for i in range(m.getnumvar())}
-            sol = pd.Series(sol, dtype=float)
+            sol_values = np.asarray(m.getxx(soltype), dtype=float)
+            if from_file:
+                sol = _solution_from_names(
+                    sol_values,
+                    [m.getvarname(i) for i in range(m.getnumvar())],
+                    self._n_vars,
+                )
+            else:
+                sol = _solution_from_labels(sol_values, self._vlabels, self._n_vars)
 
             try:
-                dual = m.gety(soltype)
-                dual = {m.getconname(i): dual[i] for i in range(m.getnumcon())}
-                dual = pd.Series(dual, dtype=float)
+                dual_values = np.asarray(m.gety(soltype), dtype=float)
+                if from_file:
+                    dual = _solution_from_names(
+                        dual_values,
+                        [m.getconname(i) for i in range(m.getnumcon())],
+                        self._n_cons,
+                    )
+                else:
+                    dual = _solution_from_labels(
+                        dual_values,
+                        self._clabels,
+                        self._n_cons,
+                    )
             except (mosek.Error, AttributeError):
                 logger.warning("Dual values of MILP couldn't be parsed")
-                dual = pd.Series(dtype=float)
+                dual = np.array([], dtype=float)
 
             return Solution(sol, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = maybe_adjust_objective_sign(solution, io_api, sense)
 
-        return Result(status, solution)
+        self.io_api = io_api
+        return self._make_result(status, solution, solver_model=m)
 
 
 class COPT(Solver[None]):
@@ -2384,57 +2735,28 @@ class COPT(Solver[None]):
         options for the given solver
     """
 
-    def __init(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
-
-    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: None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        msg = "Direct API not implemented for COPT"
-        raise NotImplementedError(msg)
+    display_name: ClassVar[str] = "COPT"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+        }
+    )
 
-    def solve_problem_from_file(
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the COPT solver.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            COPT environment for the solver
-
-        Returns
-        -------
-        Result
-        """
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         # conditions: https://guide.coap.online/copt/en-doc/constant.html#chapconst-solstatus
         CONDITION_MAP = {
             0: "unstarted",
@@ -2493,13 +2815,23 @@ def get_solver_solution() -> Solution:
             # TODO: check if this suffices
             objective = m.BestObj if m.ismip else m.LpObjVal
 
-            sol = pd.Series({v.name: v.x for v in m.getVars()}, dtype=float)
+            vars_ = m.getVars()
+            sol = _solution_from_names(
+                np.array([v.x for v in vars_], dtype=float),
+                [v.name for v in vars_],
+                self._n_vars,
+            )
 
             try:
-                dual = pd.Series({v.name: v.pi for v in m.getConstrs()}, dtype=float)
+                cons = m.getConstrs()
+                dual = _solution_from_names(
+                    np.array([c.pi for c in cons], dtype=float),
+                    [c.name for c in cons],
+                    self._n_cons,
+                )
             except (coptpy.CoptError, AttributeError):
                 logger.warning("Dual values of MILP couldn't be parsed")
-                dual = pd.Series(dtype=float)
+                dual = np.array([], dtype=float)
 
             return Solution(sol, dual, objective)
 
@@ -2508,7 +2840,8 @@ def get_solver_solution() -> Solution:
 
         env_.close()
 
-        return Result(status, solution, m)
+        self.io_api = io_api
+        return self._make_result(status, solution, solver_model=m)
 
 
 class MindOpt(Solver[None]):
@@ -2526,58 +2859,28 @@ class MindOpt(Solver[None]):
         options for the given solver
     """
 
-    def __init(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
-
-    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: None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
-    ) -> Result:
-        msg = "Direct API not implemented for MindOpt"
-        raise NotImplementedError(msg)
+    display_name: ClassVar[str] = "MindOpt"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.INTEGER_VARIABLES,
+            SolverFeature.QUADRATIC_OBJECTIVE,
+            SolverFeature.LP_FILE_NAMES,
+            SolverFeature.READ_MODEL_FROM_FILE,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+        }
+    )
 
-    def solve_problem_from_file(
+    def _run_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: None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the MindOpt solver.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            MindOpt environment for the solver
-
-        Returns
-        -------
-        Result
-
-        """
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         CONDITION_MAP = {
             -1: "error",
             0: "unknown",
@@ -2637,13 +2940,23 @@ def solve_problem_from_file(
         def get_solver_solution() -> Solution:
             objective = m.objval
 
-            sol = pd.Series({v.varname: v.X for v in m.getVars()}, dtype=float)
+            vars_ = m.getVars()
+            sol = _solution_from_names(
+                np.array([v.X for v in vars_], dtype=float),
+                [v.VarName for v in vars_],
+                self._n_vars,
+            )
 
             try:
-                dual = pd.Series({c.constrname: c.DualSoln for c in m.getConstrs()})
+                cons = m.getConstrs()
+                dual = _solution_from_names(
+                    np.array([c.DualSoln for c in cons], dtype=float),
+                    [c.ConstrName for c in cons],
+                    self._n_cons,
+                )
             except (mindoptpy.MindoptError, AttributeError):
                 logger.warning("Dual values of MILP couldn't be parsed")
-                dual = pd.Series(dtype=float)
+                dual = np.array([], dtype=float)
 
             return Solution(sol, dual, objective)
 
@@ -2653,7 +2966,8 @@ def get_solver_solution() -> Solution:
         m.dispose()
         env_.dispose()
 
-        return Result(status, solution, m)
+        self.io_api = io_api
+        return self._make_result(status, solution, solver_model=m)
 
 
 class PIPS(Solver[None]):
@@ -2661,11 +2975,7 @@ class PIPS(Solver[None]):
     Solver subclass for the PIPS solver.
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
+    def __post_init__(self) -> None:
         msg = "The PIPS solver interface is not yet implemented."
         raise NotImplementedError(msg)
 
@@ -2690,48 +3000,26 @@ class cuPDLPx(Solver[None]):
         options for the given solver
     """
 
-    def __init__(
-        self,
-        **solver_options: Any,
-    ) -> None:
-        super().__init__(**solver_options)
+    display_name: ClassVar[str] = "cuPDLPx"
+    features: ClassVar[frozenset[SolverFeature]] = frozenset(
+        {
+            SolverFeature.DIRECT_API,
+            SolverFeature.GPU_ACCELERATION,
+            SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+        }
+    )
 
-    def solve_problem_from_file(
+    def _run_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: EnvType | None = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem from a problem file using the solver cuPDLPx.
-        cuPDLPx does not currently support its own file IO, so this function
-        reads the problem file using linopy (only support netcf files) and
-        then passes the model to cuPDLPx for solving.
-        If the solution is feasible the function returns the
-        objective, solution and dual constraint variables.
-
-        Parameters
-        ----------
-        problem_fn : Path
-            Path to the problem file.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
-
-        Returns
-        -------
-        Result
-        """
+        problem_fn = self._problem_fn
+        assert problem_fn is not None
         logger.warning(
             "cuPDLPx doesn't currently support file IO. Building model from file using linopy."
         )
@@ -2743,8 +3031,9 @@ def solve_problem_from_file(
             msg = "linopy currently only supports reading models from netcdf files. Try using io_api='direct' instead."
             raise NotImplementedError(msg)
 
-        return self.solve_problem_from_model(
-            model,
+        self.model = model
+        self._build_direct()
+        return self._run_direct(
             solution_fn=solution_fn,
             log_fn=log_fn,
             warmstart_fn=warmstart_fn,
@@ -2752,67 +3041,85 @@ def solve_problem_from_file(
             env=env,
         )
 
-    def solve_problem_from_model(
+    def _build_direct(self, **kwargs: Any) -> None:
+        model = self.model
+        assert model is not None
+        if model.type in ["QP", "MILP"]:
+            msg = "cuPDLPx does not currently support QP or MILP problems."
+            raise NotImplementedError(msg)
+        if kwargs.get("explicit_coordinate_names"):
+            warnings.warn(
+                "cuPDLPx does not support named variables/constraints. "
+                "The explicit_coordinate_names parameter is ignored.",
+                UserWarning,
+                stacklevel=2,
+            )
+        cu_model = self._build_solver_model(model)
+        self.solver_model = cu_model
+        self.io_api = "direct"
+        self.sense = model.sense
+        self._cache_model_labels(model)
+
+    @staticmethod
+    def _build_solver_model(model: Model) -> cupdlpx.Model:
+        """Build a cupdlpx.Model that mirrors the linopy `model`."""
+        if model.variables.semi_continuous:
+            raise NotImplementedError(
+                "Semi-continuous variables are not supported by cuPDLPx. "
+                "Use a solver that supports them (gurobi, cplex, highs)."
+            )
+
+        M = model.matrices
+        if M.A is None:
+            raise ValueError("Model has no constraints, cannot export to cuPDLPx.")
+        A = M.A.tocsr()
+        lower = np.where(
+            np.logical_or(np.equal(M.sense, ">"), np.equal(M.sense, "=")),
+            M.b,
+            -np.inf,
+        )
+        upper = np.where(
+            np.logical_or(np.equal(M.sense, "<"), np.equal(M.sense, "=")),
+            M.b,
+            np.inf,
+        )
+
+        cu_model = cupdlpx.Model(
+            objective_vector=M.c,
+            constraint_matrix=A,
+            constraint_lower_bound=lower,
+            constraint_upper_bound=upper,
+            variable_lower_bound=M.lb,
+            variable_upper_bound=M.ub,
+        )
+
+        if model.objective.sense == "max":
+            cu_model.ModelSense = cupdlpx.PDLP.MAXIMIZE
+
+        return cu_model
+
+    def _run_direct(
         self,
-        model: Model,
         solution_fn: Path | None = None,
         log_fn: Path | None = None,
         warmstart_fn: Path | None = None,
         basis_fn: Path | None = None,
-        env: EnvType | None = None,
-        explicit_coordinate_names: bool = False,
-        set_names: bool = True,
+        env: Any = None,
+        **kw: Any,
     ) -> Result:
-        """
-        Solve a linear problem directly from a linopy model using the solver cuPDLPx.
-        If the solution is feasible the function returns the
-        objective, solution and dual constraint variables.
-
-        Parameters
-        ----------
-        model : linopy.model
-            Linopy model for the problem.
-        solution_fn : Path, optional
-            Path to the solution file.
-        log_fn : Path, optional
-            Path to the log file.
-        warmstart_fn : Path, optional
-            Path to the warmstart file.
-        basis_fn : Path, optional
-            Path to the basis file.
-        env : None, optional
-            Environment for the solver
-        explicit_coordinate_names : bool, optional
-            Transfer variable and constraint names to the solver (default: False)
-        set_names : bool, optional
-            Ignored. cuPDLPx does not support named variables/constraints.
-
-        Returns
-        -------
-        Result
-        """
-
-        if model.type in ["QP", "MILP"]:
-            msg = "cuPDLPx does not currently support QP or MILP problems."
-            raise NotImplementedError(msg)
-
-        cu_model = model.to_cupdlpx()
-
         return self._solve(
-            cu_model,
-            l_model=model,
+            self.solver_model,
             solution_fn=solution_fn,
             log_fn=log_fn,
             warmstart_fn=warmstart_fn,
             basis_fn=basis_fn,
-            io_api="direct",
-            sense=model.sense,
+            io_api=self.io_api,
+            sense=self.sense,
         )
 
     def _solve(
         self,
         cu_model: cupdlpx.Model,
-        l_model: Model | None = None,
         solution_fn: Path | None = None,
         log_fn: Path | None = None,
         warmstart_fn: Path | None = None,
@@ -2889,23 +3196,31 @@ def _solve(
 
         def get_solver_solution() -> Solution:
             objective = cu_model.ObjVal
-
-            vlabels = None if l_model is None else l_model.matrices.vlabels
-            clabels = None if l_model is None else l_model.matrices.clabels
-
-            sol = pd.Series(cu_model.X, vlabels, dtype=float)
-            dual = pd.Series(cu_model.Pi, clabels, dtype=float)
+            sol = np.asarray(cu_model.X, dtype=float)
+            dual = np.asarray(cu_model.Pi, dtype=float)
 
             if cu_model.ModelSense == cupdlpx.PDLP.MAXIMIZE:
-                dual *= -1  # flip sign of duals for max problems
+                dual = -dual
+
+            sol = _solution_from_labels(sol, self._vlabels, self._n_vars)
+            dual = _solution_from_labels(dual, self._clabels, self._n_cons)
 
             return Solution(sol, dual, objective)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = maybe_adjust_objective_sign(solution, io_api, sense)
 
-        # see https://github.com/MIT-Lu-Lab/cuPDLPx/tree/main/python#solution-attributes
-        return Result(status, solution, cu_model)
+        runtime: float | None = None
+        with contextlib.suppress(Exception):
+            runtime = float(cu_model.Runtime)
+
+        self.io_api = io_api
+        return self._make_result(
+            status,
+            solution,
+            solver_model=cu_model,
+            report=SolverReport(runtime=runtime),
+        )
 
     def _set_solver_params(self, cu_model: cupdlpx.Model) -> None:
         """
@@ -2916,3 +3231,25 @@ def _set_solver_params(self, cu_model: cupdlpx.Model) -> None:
         """
         for k, v in self.solver_options.items():
             cu_model.setParam(k, v)
+
+
+def _solver_class_for(name: str) -> type[Solver] | None:
+    try:
+        return globals().get(SolverName(name).name)
+    except ValueError:
+        return None
+
+
+QUADRATIC_SOLVERS = [
+    n.value
+    for n in SolverName
+    if (cls := _solver_class_for(n.value)) is not None
+    and cls.supports(SolverFeature.QUADRATIC_OBJECTIVE)
+]
+NO_SOLUTION_FILE_SOLVERS = [
+    n.value
+    for n in SolverName
+    if (cls := _solver_class_for(n.value)) is not None
+    and cls.supports(SolverFeature.SOLUTION_FILE_NOT_NEEDED)
+]
+quadratic_solvers = [s for s in QUADRATIC_SOLVERS if s in available_solvers]
diff --git a/linopy/variables.py b/linopy/variables.py
index dfc49a8f..cbf2fb87 100644
--- a/linopy/variables.py
+++ b/linopy/variables.py
@@ -62,7 +62,6 @@
     SOS_TYPE_ATTR,
     TERM_DIM,
 )
-from linopy.solver_capabilities import SolverFeature, solver_supports
 from linopy.types import (
     ConstantLike,
     DimsLike,
@@ -977,9 +976,11 @@ def get_solver_attribute(self, attr: str) -> DataArray:
         -------
         xr.DataArray
         """
+        from linopy.solver_capabilities import SolverFeature, solver_supports
+
         solver_model = self.model.solver_model
         if not solver_supports(
-            self.model.solver_name, SolverFeature.SOLVER_ATTRIBUTE_ACCESS
+            self.model.solver_name or "", SolverFeature.SOLVER_ATTRIBUTE_ACCESS
         ):
             raise NotImplementedError(
                 "Solver attribute getter only supports the Gurobi solver for now."
diff --git a/pyproject.toml b/pyproject.toml
index cfbaa10a..f5fa135a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -75,6 +75,7 @@ dev = [
     "types-requests",
     "gurobipy",
     "highspy",
+    "jupyter",
 ]
 benchmarks = [
     "pytest-benchmark",
diff --git a/test/test_constraint_label_index.py b/test/test_constraint_label_index.py
new file mode 100644
index 00000000..bcb506db
--- /dev/null
+++ b/test/test_constraint_label_index.py
@@ -0,0 +1,84 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+import linopy
+import linopy.constants
+from linopy.constraints import Constraint
+
+
+@pytest.fixture
+def model_with_mask() -> linopy.Model:
+    m = linopy.Model()
+    coords = pd.Index(range(5), name="i")
+    mask = np.array([True, False, True, True, False])
+    x = m.add_variables(lower=0, coords=[coords], name="x")
+    y = m.add_variables(lower=0, coords=[coords], name="y")
+    m.add_constraints(x + y >= 1, name="c_xy", mask=mask)
+    m.add_constraints(x.sum() + y.sum() <= 100, name="c_sum")
+    m.add_objective(x.sum() + 2 * y.sum())
+    return m
+
+
+def test_clabels_parity_with_matrices(model_with_mask: linopy.Model) -> None:
+    expected = model_with_mask.matrices.clabels
+    actual = model_with_mask.constraints.label_index.clabels
+    np.testing.assert_array_equal(actual, expected)
+
+
+def test_apply_result_does_not_build_matrix(
+    monkeypatch: pytest.MonkeyPatch, model_with_mask: linopy.Model
+) -> None:
+    calls = {"n": 0}
+    original = Constraint._matrix_export_data
+
+    def counting(self, label_index):  # type: ignore[no-untyped-def]
+        calls["n"] += 1
+        return original(self, label_index)
+
+    monkeypatch.setattr(Constraint, "_matrix_export_data", counting)
+
+    model_with_mask.solve("highs")
+
+    assert model_with_mask.status == "ok"
+    # one build for solver input is fine; the post-solve mapping must not add more
+    n_after_solve = calls["n"]
+    solver = model_with_mask.solver
+    result = linopy.constants.Result(
+        status=solver.status,
+        solution=solver.solution,
+        solver_model=solver.solver_model,
+        solver_name=solver.solver_name.value,
+        report=solver.report,
+    )
+    model_with_mask.apply_result(result)
+    assert calls["n"] == n_after_solve
+
+
+def test_label_index_invalidated_on_add(model_with_mask: linopy.Model) -> None:
+    first = model_with_mask.constraints.label_index.clabels.copy()
+    x = model_with_mask.variables["x"]
+    model_with_mask.add_constraints(x.sum() >= 0, name="c_extra")
+    second = model_with_mask.constraints.label_index.clabels
+    assert len(second) == len(first) + 1
+
+
+def test_label_index_invalidated_on_remove(model_with_mask: linopy.Model) -> None:
+    before = len(model_with_mask.constraints.label_index.clabels)
+    removed = len(model_with_mask.constraints["c_sum"].active_labels())
+    model_with_mask.constraints.remove("c_sum")
+    after = len(model_with_mask.constraints.label_index.clabels)
+    assert after == before - removed
+
+
+def test_apply_result_correctness_with_mask(model_with_mask: linopy.Model) -> None:
+    model_with_mask.solve("highs")
+    assert model_with_mask.status == "ok"
+    x_sol = model_with_mask.variables["x"].solution.values
+    y_sol = model_with_mask.variables["y"].solution.values
+    assert np.isfinite(x_sol).all()
+    assert np.isfinite(y_sol).all()
+    dual = model_with_mask.constraints["c_xy"].dual.values
+    mask = np.array([True, False, True, True, False])
+    assert np.isfinite(dual[mask]).all()
+    assert np.isnan(dual[~mask]).all()
diff --git a/test/test_infeasibility.py b/test/test_infeasibility.py
index feb7782d..78073acf 100644
--- a/test/test_infeasibility.py
+++ b/test/test_infeasibility.py
@@ -166,9 +166,7 @@ def test_no_solver_model_error(self, solver: str) -> None:
         # Solve the model first
         m.solve(solver_name=solver)
 
-        # Manually remove the solver_model to simulate cleanup
-        m.solver_model = None
-        m.solver_name = solver  # But keep the solver name
+        m.solver.solver_model = None
 
         # Should raise ValueError since we know it was solved with supported solver
         with pytest.raises(ValueError, match="No solver model available"):
diff --git a/test/test_io.py b/test/test_io.py
index 0a4c4e64..b049c0dc 100644
--- a/test/test_io.py
+++ b/test/test_io.py
@@ -274,7 +274,8 @@ def test_to_file_invalid(model: Model, tmp_path: Path) -> None:
 
 @pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobipy not installed")
 def test_to_gurobipy(model: Model) -> None:
-    model.to_gurobipy()
+    gm = model.to_gurobipy()
+    assert gm.NumVars > 0
 
 
 @pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobipy not installed")
@@ -288,7 +289,8 @@ def test_to_gurobipy_no_names(model: Model) -> None:
 
 @pytest.mark.skipif("highs" not in available_solvers, reason="Highspy not installed")
 def test_to_highspy(model: Model) -> None:
-    model.to_highspy()
+    h = model.to_highspy()
+    assert h.getLp().num_col_ > 0
 
 
 @pytest.mark.skipif("highs" not in available_solvers, reason="Highspy not installed")
@@ -299,6 +301,18 @@ def test_to_highspy_no_names(model: Model) -> None:
     assert len(lp.row_names_) == 0
 
 
+@pytest.mark.skipif("mosek" not in available_solvers, reason="Mosek not installed")
+def test_to_mosek(model: Model) -> None:
+    task = model.to_mosek()
+    assert task.getnumvar() > 0
+
+
+@pytest.mark.skipif("cupdlpx" not in available_solvers, reason="cuPDLPx not installed")
+def test_to_cupdlpx(model: Model) -> None:
+    cu = model.to_cupdlpx()
+    assert cu is not None
+
+
 def test_model_set_names_in_solver_io_default() -> None:
     assert Model().set_names_in_solver_io is True
 
diff --git a/test/test_optimization.py b/test/test_optimization.py
index 07d23fa2..9b3eeb09 100644
--- a/test/test_optimization.py
+++ b/test/test_optimization.py
@@ -25,7 +25,13 @@
     get_available_solvers_with_feature,
     solver_supports,
 )
-from linopy.solvers import _new_highspy_mps_layout, available_solvers, quadratic_solvers
+from linopy.solvers import (
+    SolverFeature,
+    _new_highspy_mps_layout,
+    _solver_class_for,
+    available_solvers,
+    quadratic_solvers,
+)
 
 logger = logging.getLogger(__name__)
 
@@ -496,6 +502,22 @@ def test_mock_solve(model_maximization: Model) -> None:
     assert (x_solution == 0).all()
 
 
+@pytest.mark.skipif("highs" not in available_solvers, reason="HiGHS is not installed")
+def test_mock_solve_clears_existing_solver_state(model: Model) -> None:
+    status, condition = model.solve(solver_name="highs", io_api="direct")
+    assert status == "ok"
+    assert model.solver is not None
+    assert model.solver_model is not None
+    assert model.solver_name == "highs"
+
+    status, condition = model.solve(solver="some_non_existant_solver", mock_solve=True)
+
+    assert status == "ok"
+    assert model.solver is None
+    assert model.solver_model is None
+    assert model.solver_name is None
+
+
 @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params)
 def test_default_settings_chunked(
     model_chunked: Model, solver: str, io_api: str, explicit_coordinate_names: bool
@@ -754,6 +776,14 @@ def test_milp_model(
     assert condition == "optimal"
     assert ((milp_model.solution.y == 9) | (milp_model.solution.x == 0.5)).all()
 
+    solver_cls = _solver_class_for(solver)
+    if solver_cls is not None and solver_cls.supports(
+        SolverFeature.MIP_DUAL_BOUND_REPORT
+    ):
+        report = milp_model.solver.report
+        assert report is not None
+        assert report.dual_bound is not None
+
 
 @pytest.mark.parametrize(
     "solver,io_api,explicit_coordinate_names",
diff --git a/test/test_solution_lookup.py b/test/test_solution_lookup.py
index 7dd9643f..3f6475a8 100644
--- a/test/test_solution_lookup.py
+++ b/test/test_solution_lookup.py
@@ -1,73 +1,54 @@
 import numpy as np
-import pandas as pd
 from numpy import nan
 
-from linopy.common import lookup_vals, series_to_lookup_array
+from linopy.common import values_to_lookup_array
+from linopy.solvers import _solution_from_names
 
 
-class TestSeriesToLookupArray:
+class TestValuesToLookupArray:
     def test_basic(self) -> None:
-        s = pd.Series([10.0, 20.0, 30.0], index=pd.Index([0, 1, 2]))
-        arr = series_to_lookup_array(s)
+        arr = values_to_lookup_array(np.array([10.0, 20.0, 30.0]), np.array([0, 1, 2]))
         np.testing.assert_array_equal(arr, [10.0, 20.0, 30.0])
 
-    def test_with_negative_index(self) -> None:
-        s = pd.Series([nan, 10.0, 20.0], index=pd.Index([-1, 0, 2]))
-        arr = series_to_lookup_array(s)
+    def test_negative_labels_skipped(self) -> None:
+        arr = values_to_lookup_array(np.array([nan, 10.0, 20.0]), np.array([-1, 0, 2]))
         assert arr[0] == 10.0
         assert np.isnan(arr[1])
         assert arr[2] == 20.0
 
-    def test_sparse_index(self) -> None:
-        s = pd.Series([5.0, 7.0], index=pd.Index([0, 100]))
-        arr = series_to_lookup_array(s)
+    def test_sparse_labels(self) -> None:
+        arr = values_to_lookup_array(np.array([5.0, 7.0]), np.array([0, 100]))
         assert len(arr) == 101
         assert arr[0] == 5.0
         assert arr[100] == 7.0
         assert np.isnan(arr[50])
 
-    def test_only_negative_index(self) -> None:
-        s = pd.Series([nan], index=pd.Index([-1]))
-        arr = series_to_lookup_array(s)
-        assert len(arr) == 1
-        assert np.isnan(arr[0])
+    def test_only_negative_labels(self) -> None:
+        arr = values_to_lookup_array(np.array([nan]), np.array([-1]))
+        assert len(arr) == 0
 
-
-class TestLookupVals:
-    def test_basic(self) -> None:
-        arr = np.array([10.0, 20.0, 30.0])
-        idx = np.array([0, 1, 2])
-        result = lookup_vals(arr, idx)
-        np.testing.assert_array_equal(result, [10.0, 20.0, 30.0])
-
-    def test_negative_labels_become_nan(self) -> None:
-        arr = np.array([10.0, 20.0])
-        idx = np.array([0, -1, 1, -1])
-        result = lookup_vals(arr, idx)
-        assert result[0] == 10.0
-        assert np.isnan(result[1])
-        assert result[2] == 20.0
-        assert np.isnan(result[3])
-
-    def test_out_of_range_labels_become_nan(self) -> None:
-        arr = np.array([10.0, 20.0])
-        idx = np.array([0, 1, 999])
-        result = lookup_vals(arr, idx)
-        assert result[0] == 10.0
-        assert result[1] == 20.0
-        assert np.isnan(result[2])
-
-    def test_all_negative(self) -> None:
-        arr = np.array([10.0])
-        idx = np.array([-1, -1, -1])
-        result = lookup_vals(arr, idx)
-        assert all(np.isnan(result))
-
-    def test_no_mutation_of_source(self) -> None:
-        arr = np.array([10.0, 20.0, 30.0])
-        idx1 = np.array([-1, 1])
-        idx2 = np.array([0, 2])
-        lookup_vals(arr, idx1)
-        result2 = lookup_vals(arr, idx2)
-        np.testing.assert_array_equal(result2, [10.0, 30.0])
-        np.testing.assert_array_equal(arr, [10.0, 20.0, 30.0])
+    def test_explicit_size(self) -> None:
+        arr = values_to_lookup_array(np.array([5.0, 7.0]), np.array([0, 2]), size=5)
+        assert len(arr) == 5
+        assert arr[0] == 5.0
+        assert arr[2] == 7.0
+        assert np.isnan(arr[1])
+        assert np.isnan(arr[3])
+        assert np.isnan(arr[4])
+
+
+class TestSolutionFromNames:
+    def test_default_names(self) -> None:
+        arr = _solution_from_names(
+            np.array([1.0, 2.0, 3.0]), ["x2", "x0", "x1"], size=4
+        )
+        np.testing.assert_array_equal(arr[:3], [2.0, 3.0, 1.0])
+        assert np.isnan(arr[3])
+
+    def test_explicit_coordinate_names(self) -> None:
+        arr = _solution_from_names(
+            np.array([1.0, 2.0]), ["power[1]#5", "power[0]#3"], size=7
+        )
+        assert arr[3] == 2.0
+        assert arr[5] == 1.0
+        assert np.isnan(arr[4])
diff --git a/test/test_solvers.py b/test/test_solvers.py
index 7f4d55ec..68e31fa4 100644
--- a/test/test_solvers.py
+++ b/test/test_solvers.py
@@ -7,11 +7,201 @@
 
 from pathlib import Path
 
+import numpy as np
 import pytest
 from test_io import model  # noqa: F401
 
-from linopy import Model, solvers
-from linopy.solver_capabilities import SolverFeature, solver_supports
+from linopy import GREATER_EQUAL, Model, solvers
+from linopy.constants import Result, Solution, Status
+from linopy.constraints import CSRConstraint
+from linopy.solver_capabilities import (
+    SOLVER_REGISTRY,
+    SolverFeature,
+    SolverInfo,
+    solver_supports,
+)
+from linopy.solvers import _installed_version_in
+
+
+@pytest.fixture
+def simple_model() -> Model:
+    m = Model(chunk=None)
+    x = m.add_variables(name="x")
+    y = m.add_variables(name="y")
+    m.add_constraints(2 * x + 6 * y, GREATER_EQUAL, 10)
+    m.add_constraints(4 * x + 2 * y, GREATER_EQUAL, 3)
+    m.add_objective(2 * y + x)
+    return m
+
+
+@pytest.mark.parametrize("solver", sorted(set(solvers.available_solvers)))
+def test_solver_instance_attached_after_solve(simple_model: Model, solver: str) -> None:
+    simple_model.solve(solver)
+    assert isinstance(simple_model.solver, solvers.Solver)
+    assert simple_model.solver.status is not None
+    assert simple_model.solver.status.is_ok
+    assert simple_model.solver.solution is not None
+    assert simple_model.solver_model is simple_model.solver.solver_model
+    assert simple_model.solver_name == solver
+
+
+@pytest.mark.parametrize("solver", sorted(set(solvers.available_solvers)))
+def test_result_carries_solver_name(simple_model: Model, solver: str) -> None:
+    if not solver_supports(solver, SolverFeature.DIRECT_API):
+        pytest.skip("Solver does not support direct API.")
+    instance = solvers.Solver.from_name(solver, simple_model, io_api="direct")
+    result = instance.solve()
+    assert result.solver_name == solver
+
+
+@pytest.mark.parametrize("solver", sorted(set(solvers.available_solvers)))
+def test_from_name_then_solve(simple_model: Model, solver: str) -> None:
+    if not solver_supports(solver, SolverFeature.DIRECT_API):
+        pytest.skip("Solver does not support direct API.")
+    built = solvers.Solver.from_name(solver, simple_model, io_api="direct")
+    assert built.solver_model is not None
+    result = built.solve()
+    simple_model.apply_result(result)
+
+    reference = Model(chunk=None)
+    rx = reference.add_variables(name="x")
+    ry = reference.add_variables(name="y")
+    reference.add_constraints(2 * rx + 6 * ry, GREATER_EQUAL, 10)
+    reference.add_constraints(4 * rx + 2 * ry, GREATER_EQUAL, 3)
+    reference.add_objective(2 * ry + rx)
+    reference.solve(solver, io_api="direct")
+
+    assert simple_model.status == "ok"
+    assert np.isclose(simple_model.objective.value, reference.objective.value)
+
+
+@pytest.mark.parametrize("solver", sorted(set(solvers.available_solvers)))
+def test_from_name_set_names_false(simple_model: Model, solver: str) -> None:
+    if not solver_supports(solver, SolverFeature.DIRECT_API):
+        pytest.skip("Solver does not support direct API.")
+    built = solvers.Solver.from_name(
+        solver, simple_model, io_api="direct", set_names=False
+    )
+    result = built.solve()
+    status, condition = simple_model.apply_result(result)
+
+    assert status == "ok"
+    assert condition == "optimal"
+    assert simple_model.objective.value == pytest.approx(3.3)
+    assert float(simple_model.variables["x"].solution) == pytest.approx(-0.1)
+    assert float(simple_model.variables["y"].solution) == pytest.approx(1.7)
+
+
+def test_from_name_unknown_solver_raises(simple_model: Model) -> None:
+    with pytest.raises(ValueError, match="unknown solver"):
+        solvers.Solver.from_name("not_a_real_solver", simple_model, io_api="direct")
+
+
+@pytest.mark.skipif(
+    "highs" not in set(solvers.available_solvers), reason="HiGHS is not installed"
+)
+def test_from_name_applies_solver_options(simple_model: Model) -> None:
+    built = solvers.Solver.from_name(
+        "highs", simple_model, io_api="direct", options={"time_limit": 123}
+    )
+    option_status, time_limit = built.solver_model.getOptionValue("time_limit")
+    assert str(option_status) == "HighsStatus.kOk"
+    assert time_limit == 123
+
+
+@pytest.mark.skipif(
+    "highs" not in set(solvers.available_solvers), reason="HiGHS is not installed"
+)
+def test_solver_state_compatibility_setters(simple_model: Model) -> None:
+    simple_model.solver = solvers.Solver.from_name(
+        "highs", simple_model, io_api="direct"
+    )
+    simple_model.solver_model = None
+    assert simple_model.solver is None
+    assert simple_model.solver_model is None
+    assert simple_model.solver_name is None
+
+    simple_model.solver = solvers.Solver.from_name(
+        "highs", simple_model, io_api="direct"
+    )
+    simple_model.solver_name = None
+    assert simple_model.solver is None
+    assert simple_model.solver_model is None
+    assert simple_model.solver_name is None
+
+    with pytest.raises(AttributeError, match="managed via model.solver"):
+        simple_model.solver_model = object()
+    with pytest.raises(AttributeError, match="managed via model.solver"):
+        simple_model.solver_name = "highs"
+
+
+def test_apply_result_explicit(simple_model: Model) -> None:
+    x_labels = simple_model.variables["x"].labels.values
+    y_labels = simple_model.variables["y"].labels.values
+    primal = np.full(simple_model._xCounter, np.nan)
+    primal[int(x_labels)] = 1.5
+    primal[int(y_labels)] = 2.0
+    solution = Solution(primal=primal, objective=5.5)
+    result = Result(
+        status=Status.from_termination_condition("optimal"),
+        solution=solution,
+        solver_name="mock",
+    )
+    simple_model.solver = None
+    simple_model.apply_result(result)
+    assert simple_model.status == "ok"
+    assert simple_model.termination_condition == "optimal"
+    assert simple_model.objective.value == 5.5
+    assert float(simple_model.variables["x"].solution) == 1.5
+    assert float(simple_model.variables["y"].solution) == 2.0
+
+
+def test_apply_result_with_csr_constraints_avoids_data_reconstruction(
+    monkeypatch: pytest.MonkeyPatch,
+) -> None:
+    m = Model(freeze_constraints=True)
+    x = m.add_variables(coords=[range(3)], name="x")
+    m.add_constraints(x >= 0, name="c")
+    con = m.constraints["c"]
+    assert isinstance(con, CSRConstraint)
+
+    primal = np.arange(m._xCounter, dtype=float)
+    dual = np.arange(m._cCounter, dtype=float) + 10
+    result = Result(
+        status=Status.from_termination_condition("optimal"),
+        solution=Solution(primal=primal, dual=dual, objective=1.0),
+        solver_name="mock",
+    )
+
+    def fail_data(self: CSRConstraint) -> None:
+        raise AssertionError("CSRConstraint.data was accessed")
+
+    monkeypatch.setattr(CSRConstraint, "data", property(fail_data))
+    m.apply_result(result)
+
+    np.testing.assert_array_equal(m.variables["x"].solution.values, primal)
+    np.testing.assert_array_equal(m.constraints["c"].dual.values, dual)
+
+
+@pytest.mark.skipif(
+    "gurobi" not in set(solvers.available_solvers), reason="Gurobi is not installed"
+)
+def test_gurobi_env_persists_after_solve(simple_model: Model) -> None:
+    simple_model.solve("gurobi", io_api="direct")
+    assert simple_model.solver is not None
+    assert simple_model.solver.env is not None
+    assert isinstance(simple_model.solver_model.NumVars, int)
+
+
+@pytest.mark.parametrize("solver", sorted(set(solvers.available_solvers)))
+def test_solver_close_releases_state(simple_model: Model, solver: str) -> None:
+    simple_model.solve(solver)
+    solver_instance = simple_model.solver
+    assert solver_instance is not None
+    solver_instance.close()
+    assert solver_instance.solver_model is None
+    assert solver_instance.env is None
+
 
 free_mps_problem = """NAME               sample_mip
 ROWS
@@ -124,7 +314,7 @@ def test_knitro_solver_for_lp(tmp_path: Path) -> None:
 )
 def test_knitro_solver_with_options(tmp_path: Path) -> None:
     """Test Knitro solver with custom options."""
-    knitro = solvers.Knitro(maxit=100, feastol=1e-6)
+    knitro = solvers.Knitro(options={"maxit": 100, "feastol": 1e-6})
 
     mps_file = tmp_path / "problem.mps"
     mps_file.write_text(free_mps_problem)
@@ -154,7 +344,7 @@ def test_knitro_solver_with_model_raises_error(model: Model) -> None:  # noqa: F
 )
 def test_knitro_solver_no_log(tmp_path: Path) -> None:
     """Test Knitro solver without log file."""
-    knitro = solvers.Knitro(outlev=0)
+    knitro = solvers.Knitro(options={"outlev": 0})
 
     mps_file = tmp_path / "problem.mps"
     mps_file.write_text(free_mps_problem)
@@ -218,3 +408,57 @@ def test_gurobi_environment_with_gurobi_env(model: Model, tmp_path: Path) -> Non
         gurobi.solve_problem(model=model, solution_fn=sol_file, env=env)
     assert result.status.is_ok
     assert log2_file.exists()
+
+
+@pytest.mark.parametrize(
+    "solver_cls, feature, expected",
+    [
+        (solvers.Gurobi, SolverFeature.SOS_CONSTRAINTS, True),
+        (solvers.Gurobi, SolverFeature.GPU_ACCELERATION, False),
+        (solvers.Highs, SolverFeature.SOS_CONSTRAINTS, False),
+        (solvers.Highs, SolverFeature.SEMI_CONTINUOUS_VARIABLES, True),
+        (solvers.CBC, SolverFeature.LP_FILE_NAMES, False),
+        (solvers.CBC, SolverFeature.INTEGER_VARIABLES, True),
+        (solvers.cuPDLPx, SolverFeature.DIRECT_API, True),
+        (solvers.cuPDLPx, SolverFeature.GPU_ACCELERATION, True),
+        (solvers.cuPDLPx, SolverFeature.QUADRATIC_OBJECTIVE, False),
+        (solvers.PIPS, SolverFeature.INTEGER_VARIABLES, False),
+    ],
+)
+def test_solver_class_supports_feature(
+    solver_cls: type, feature: SolverFeature, expected: bool
+) -> None:
+    assert solver_cls.supports(feature) is expected
+
+
+def test_solver_instance_supports_matches_class() -> None:
+    feature = SolverFeature.QUADRATIC_OBJECTIVE
+    assert solvers.Gurobi.supports(feature) is True
+    if "gurobi" in solvers.available_solvers:
+        assert solvers.Gurobi().supports(feature) is True
+
+
+@pytest.mark.parametrize("solver_name", [n.value for n in solvers.SolverName])
+def test_capability_shim_round_trips(solver_name: str) -> None:
+    solver_cls = getattr(solvers, solvers.SolverName(solver_name).name)
+    for feature in SolverFeature:
+        assert solver_supports(solver_name, feature) == solver_cls.supports(feature)
+
+
+def test_solver_registry_iter_and_index() -> None:
+    names = list(SOLVER_REGISTRY)
+    assert "gurobi" in names
+    for name in names:
+        info = SOLVER_REGISTRY[name]
+        assert isinstance(info, SolverInfo)
+        assert isinstance(info.features, frozenset)
+        assert info.name == name
+
+
+@pytest.mark.skipif(
+    "xpress" not in set(solvers.available_solvers), reason="Xpress is not installed"
+)
+def test_xpress_gpu_feature_reflects_installed_version() -> None:
+    assert solvers.Xpress.supports(
+        SolverFeature.GPU_ACCELERATION
+    ) == _installed_version_in("xpress", ">=9.8.0")