Improve Mosek interface to select optimal solution

doc/release_notes.rst

@@ -54,6 +54,7 @@ Most users should keep calling ``model.solve(...)``. If you want more control, y
 
 * SOS constraints on masked variables no longer cause solver-specific failures (Gurobi ``IndexError``, Xpress ``?404 Invalid column number``, LP parse errors, silent set corruption). ``Model.solve()`` and ``Model.to_file()`` now raise a clear ``NotImplementedError`` referring users to `#688 <https://github.com/PyPSA/linopy/issues/688>`__; pass ``reformulate_sos=True`` as a workaround.
 * ``Model.solve(..., reformulate_sos=True)`` now actually reformulates SOS constraints even when the solver supports them natively. Previously it was silently ignored with a warning.
+* Fix Mosek interface to inspect both the basic and IPM solutions and pick the one with the better status, so that an optimal crossover solution is not discarded when IPM terminates with a (near-)Farkas certificate.
 
 **Breaking Changes**
 

linopy/solvers.py

@@ -2868,6 +2868,58 @@ def _build_solver_model(
             task.putobjsense(mosek.objsense.minimize)
         return task
 
+    @staticmethod
+    def _choose_solution(task: mosek.Task) -> mosek.soltype | None:
+        """
+        Pick the Mosek solution with the best status available.
+
+        Mosek may return up to three solutions per task: interior-point
+        (``soltype.itr``), basic (``soltype.bas``), and integer
+        (``soltype.itg``). Each carries its own ``solsta``: on a numerically
+        marginal LP solved with the default IPM+crossover, the interior-point
+        solver may terminate with ``solsta.dual_infeas_cer`` while crossover
+        recovers ``solsta.optimal`` for the basic solution. Reading only the
+        interior-point solution would discard the actual optimum.
+
+        Ranking, best to worst: ``solsta.optimal`` / ``solsta.integer_optimal``
+        > any other defined status > undefined. On a tie between ``bas`` and
+        ``itr`` (e.g. both ``optimal``) we prefer ``itr`` to preserve historical
+        behaviour. If ``itg`` is defined it always wins, since integer and
+        continuous solutions do not coexist for a well-posed task.
+
+        Returns ``None`` if no solution is defined at all (e.g. the optimizer
+        crashed before producing one).
+        """
+
+        def _is_defined(soltype: mosek.soltype) -> bool:
+            try:
+                return bool(task.solutiondef(soltype))
+            except mosek.Error:
+                return False
+
+        if _is_defined(mosek.soltype.itg):
+            return mosek.soltype.itg
+
+        optimal_statuses = {mosek.solsta.optimal, mosek.solsta.integer_optimal}
+
+        best: mosek.soltype | None = None
+        best_score = -1
+        # Iterate bas first and only then itr so that on a score tie
+        # itr wins, preserving the historical default for the common LP case.
+        for candidate in [mosek.soltype.bas, mosek.soltype.itr]:
+            if not _is_defined(candidate):
+                continue
+            try:
+                solsta = task.getsolsta(candidate)
+            except mosek.Error:
+                continue
+            score = 1 if solsta in optimal_statuses else 0
+            if score >= best_score:
+                best = candidate
+                best_score = score
+
+        return best
+
     def _run_file(
         self,
         solution_fn: Path | None = None,
@@ -3052,25 +3104,25 @@ def _solve(
                             f.write(f" UL {namex}\n")
                     f.write("ENDATA\n")
 
-        soltype = None
-        possible_soltypes = [
-            mosek.soltype.bas,
-            mosek.soltype.itr,
-            mosek.soltype.itg,
-        ]
-        for possible_soltype in possible_soltypes:
-            try:
-                if m.solutiondef(possible_soltype):
-                    soltype = possible_soltype
-            except mosek.Error:
-                pass
+        # Inspect both bas and itr (and itg for MILPs) and pick the
+        # solution with the best status. Reading only the interior-point
+        # solution may discard a valid crossover optimum.
+        soltype = Mosek._choose_solution(m)
 
-        if solution_fn is not None:
+        if solution_fn is not None and soltype is not None:
             try:
-                m.writesolution(mosek.soltype.bas, path_to_string(solution_fn))
+                m.writesolution(soltype, path_to_string(solution_fn))
             except mosek.Error as err:
                 logger.info("Unable to save solution file. Raised error: %s", err)
 
+        if soltype is None:
+            condition = "no solution available"
+            status = Status.from_termination_condition(
+                TerminationCondition.internal_solver_error
+            )
+            status.legacy_status = condition
+            return self._make_result(status, None)
+
         condition = str(m.getsolsta(soltype))
         termination_condition = CONDITION_MAP.get(condition, condition)
         status = Status.from_termination_condition(termination_condition)

test/test_solvers.py

@@ -6,6 +6,7 @@
 """
 
 from pathlib import Path
+from unittest.mock import MagicMock
 
 import numpy as np
 import pytest
@@ -567,3 +568,124 @@ def test_assign_result_without_solver_kwarg_leaves_solver_unset(self) -> None:
         m.assign_result(result)  # no solver kwarg
 
         assert m.solver is None
+
+
+mosek_installed = pytest.importorskip("mosek", reason="Mosek is not installed")
+
+
+class TestMosekChooseSolution:
+    @staticmethod
+    def _make_task_mock(
+        *,
+        bas_solsta: object | None = None,
+        itr_solsta: object | None = None,
+        itg_solsta: object | None = None,
+    ) -> MagicMock:
+        defined = {
+            mosek_installed.soltype.bas: bas_solsta,
+            mosek_installed.soltype.itr: itr_solsta,
+            mosek_installed.soltype.itg: itg_solsta,
+        }
+        task = MagicMock()
+        task.solutiondef.side_effect = lambda st: defined[st] is not None
+        task.getsolsta.side_effect = lambda st: defined[st]
+        return task
+
+    @pytest.mark.parametrize(
+        "kwargs, expected_soltype",
+        [
+            pytest.param(
+                dict(
+                    bas_solsta=mosek_installed.solsta.optimal,
+                    itr_solsta=mosek_installed.solsta.dual_infeas_cer,
+                ),
+                mosek_installed.soltype.bas,
+                id="prefers_bas_when_itr_is_farkas",
+            ),
+            pytest.param(
+                dict(
+                    bas_solsta=mosek_installed.solsta.optimal,
+                    itr_solsta=mosek_installed.solsta.optimal,
+                ),
+                mosek_installed.soltype.itr,
+                id="prefers_itr_on_tie",
+            ),
+            pytest.param(
+                dict(itr_solsta=mosek_installed.solsta.optimal),
+                mosek_installed.soltype.itr,
+                id="only_itr_defined",
+            ),
+            pytest.param(
+                dict(bas_solsta=mosek_installed.solsta.optimal),
+                mosek_installed.soltype.bas,
+                id="only_bas_defined",
+            ),
+            pytest.param(
+                dict(),
+                None,
+                id="nothing_defined",
+            ),
+            pytest.param(
+                dict(itg_solsta=mosek_installed.solsta.integer_optimal),
+                mosek_installed.soltype.itg,
+                id="itg_for_mip",
+            ),
+            pytest.param(
+                dict(
+                    bas_solsta=mosek_installed.solsta.optimal,
+                    itr_solsta=mosek_installed.solsta.optimal,
+                    itg_solsta=mosek_installed.solsta.integer_optimal,
+                ),
+                mosek_installed.soltype.itg,
+                id="itg_wins_over_bas_itr",
+            ),
+            pytest.param(
+                dict(
+                    bas_solsta=mosek_installed.solsta.unknown,
+                    itr_solsta=mosek_installed.solsta.optimal,
+                ),
+                mosek_installed.soltype.itr,
+                id="optimal_itr_over_unknown_bas",
+            ),
+            pytest.param(
+                dict(
+                    bas_solsta=mosek_installed.solsta.optimal,
+                    itr_solsta=mosek_installed.solsta.unknown,
+                ),
+                mosek_installed.soltype.bas,
+                id="optimal_bas_over_unknown_itr",
+            ),
+            pytest.param(
+                dict(
+                    bas_solsta=mosek_installed.solsta.prim_infeas_cer,
+                    itr_solsta=mosek_installed.solsta.dual_infeas_cer,
+                ),
+                mosek_installed.soltype.itr,
+                id="falls_back_to_itr_when_both_non_optimal",
+            ),
+        ],
+    )
+    def test_choose_solution(
+        self, kwargs: dict[str, object], expected_soltype: object
+    ) -> None:
+        task = self._make_task_mock(**kwargs)
+        assert solvers.Mosek._choose_solution(task) is expected_soltype
+
+    @pytest.mark.skipif(
+        "mosek" not in set(solvers.licensed_solvers),
+        reason="Mosek is not licensed",
+    )
+    def test_smoke_lp(self) -> None:
+        import math
+
+        m = Model()
+        x = m.add_variables(name="x", lower=0)
+        m.add_constraints(2 * x >= 10, name="c1")
+        m.add_objective(x)
+
+        result = solvers.Solver.from_name("mosek", m).solve()
+
+        assert result.status.is_ok
+        assert result.solution is not None
+        assert math.isfinite(result.solution.objective)
+        assert result.solution.objective == pytest.approx(5.0, abs=1e-3)