Add quadratic objective support and 'direct' mode for MOSEK

linopy/io.py

@@ -13,7 +13,7 @@
 import pandas as pd
 import xarray as xr
 from numpy import asarray, concatenate, ones_like, zeros_like
-from scipy.sparse import csc_matrix, triu
+from scipy.sparse import csc_matrix, tril, triu
 from tqdm import tqdm
 
 from linopy import solvers
@@ -307,6 +307,93 @@
     return fn
 
 
+def to_mosekpy(model, task=None):
+    """
+    Export model to MOSEK.
+
+    Export the model directly to MOSEK without writing files.
+
+    Parameters
+    ----------
+    m : linopy.Model
+    task : empty MOSEK task
+
+    Returns
+    -------
+    task : MOSEK Task object
+    """
+
+    import mosek
+
+    if task is None:
+        task = mosek.Task()
+
+    task.appendvars(model.nvars)
+    task.appendcons(model.ncons)
+
+    M = model.matrices
+    for j, n in enumerate(("x" + M.vlabels.astype(str).astype(object))):
+        task.putvarname(j, n)
+
+    ## Variables
+
+    bkx = [
+        (
+            (
+                (mosek.boundkey.ra if l < u else mosek.boundkey.fx)
+                if u < np.inf
+                else mosek.boundkey.lo
+            )
+            if (l > -np.inf)
+            else (mosek.boundkey.up if (u < np.inf) else mosek.boundkey.fr)
+        )
+        for (l, u) 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)
+
+    ## Constraints
+
+    if len(model.constraints) > 0:
+        names = "c" + M.clabels.astype(str).astype(object)
+        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
+        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)
+
+    ## Objective
+    if model.is_quadratic:
+        Q = (0.5 * tril(M.Q + M.Q.transpose())).tocoo()
+        task.putqobj(Q.row, Q.col, Q.data)
+    task.putclist(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 to_gurobipy(m, env=None):
     """
     Export the model to gurobipy.

linopy/model.py

@@ -27,7 +27,14 @@
     QuadraticExpression,
     ScalarLinearExpression,
 )
-from linopy.io import to_block_files, to_file, to_gurobipy, to_highspy, to_netcdf
+from linopy.io import (
+    to_block_files,
+    to_file,
+    to_gurobipy,
+    to_highspy,
+    to_mosekpy,
+    to_netcdf,
+)
 from linopy.matrices import MatrixAccessor
 from linopy.objective import Objective
 from linopy.solvers import available_solvers, quadratic_solvers
@@ -1185,6 +1192,8 @@ def reset_solution(self):
 
     to_gurobipy = to_gurobipy
 
+    to_mosekpy = to_mosekpy
+
     to_highspy = to_highspy
 
     to_block_files = to_block_files

linopy/solvers.py

@@ -31,6 +31,7 @@
     "highs",
     "scip",
     "mosek",
+    "mosek_remote",
     "copt",
     "mindopt",
 ]
@@ -77,6 +78,15 @@
             m.checkinall()
 
         available_solvers.append("mosek")
+with contextlib.suppress(ImportError):
+    import mosek
+
+    with contextlib.suppress(mosek.Error):
+        with mosek.Task() as t:
+            t.putoptserverhost("http://solve.mosek.com:30080")
+            t.optimize()
+        available_solvers.append("mosek_remote")
+
 with contextlib.suppress(ImportError):
     import mindoptpy
 
@@ -94,7 +104,15 @@
 
 
 io_structure = dict(
-    lp_file={"gurobi", "xpress", "cbc", "glpk", "cplex", "mosek", "mindopt"},
+    lp_file={
+        "gurobi",
+        "xpress",
+        "cbc",
+        "glpk",
+        "cplex",
+        "mosek",
+        "mosek_remote" "mindopt",
+    },
     blocks={"pips"},
 )
 
@@ -845,6 +863,45 @@
     return Result(status, solution, m)
 
 
+mosek_bas_re = re.compile(r" (XL|XU)\s+([^ \t]+)\s+([^ \t]+)| (LL|UL|BS)\s+([^ \t]+)")
+
+
+def run_mosek_remote(
+    model,
+    io_api=None,
+    problem_fn=None,
+    solution_fn=None,
+    log_fn=None,
+    warmstart_fn=None,
+    basis_fn=None,
+    keep_files=False,
+    env=None,
+    **solver_options,
+):
+    """
+    This is a wrapper around run_mosek() that sets the
+    MSK_SPAR_REMOTE_OPTSERVER_HOST, if not already present, to the public MOSEK
+    solver server. This can be used to solve small-ish problems without a
+    license.
+    """
+    if "MSK_SPAR_REMOTE_OPTSERVER_HOST" not in solver_options:
+        solver_options["MSK_SPAR_REMOTE_OPTSERVER_HOST"] = (
+            "http://solve.mosek.com:30080"
+        )
+    return run_mosek(
+        model,
+        io_api=io_api,
+        problem_fn=problem_fn,
+        solution_fn=solution_fn,
+        log_fn=log_fn,
+        warmstart_fn=warmstart_fn,
+        basis_fn=basis_fn,
+        keep_files=keep_files,
+        env=env,
+        **solver_options,
+    )
+
+
 def run_mosek(
     model,
     io_api=None,
@@ -858,7 +915,9 @@
     **solver_options,
 ):
     """
-    Solve a linear problem using the MOSEK solver.
+    Solve a linear problem using the MOSEK solver. Both 'direct' mode, mps and
+    lp mode are supported; None is interpret as 'direct' mode. MPS mode does
+    not support quadratic terms.
 
     https://www.mosek.com/
 
@@ -873,24 +932,28 @@
         "solsta.dual_infeas_cer": "infeasible_or_unbounded",
     }
 
-    if io_api is not None and io_api not in ["lp", "mps"]:
+    if io_api is not None and io_api not in ["lp", "mps", "direct"]:
         raise ValueError(
             "Keyword argument `io_api` has to be one of `lp`, `mps` or None"
         )
 
     problem_fn = model.to_file(problem_fn)
 
-    problem_fn = maybe_convert_path(problem_fn)
+    if io_api != "direct" and io_api is not None:
+        problem_fn = maybe_convert_path(problem_fn)
     log_fn = maybe_convert_path(log_fn)
-    warmstart_fn = maybe_convert_path(warmstart_fn)
-    basis_fn = maybe_convert_path(basis_fn)
+    # warmstart_fn = maybe_convert_path(warmstart_fn)
+    # basis_fn = maybe_convert_path(basis_fn)
 
     with contextlib.ExitStack() as stack:
         if env is None:
             env = stack.enter_context(mosek.Env())
 
         with env.Task() as m:
-            m.readdata(problem_fn)
+            if io_api is None or io_api == "direct":
+                model.to_mosekpy(m)
+            else:
+                m.readdata(problem_fn)
 
             for k, v in solver_options.items():
                 m.putparam(k, str(v))
@@ -899,17 +962,107 @@
                 m.linkfiletostream(mosek.streamtype.log, log_fn, 0)
 
             if warmstart_fn:
-                m.readdata(warmstart_fn)
-
+                m.putintparam(mosek.iparam.sim_hotstart, mosek.simhotstart.status_keys)
+                skx = [mosek.stakey.low] * m.getnumvar()
+                skc = [mosek.stakey.bas] * m.getnumcon()
+
+                with open(warmstart_fn, "rt") as f:
+                    for line in f:
+                        if line.startswith("NAME "):
+                            break
+
+                    for line in f:
+                        if line.startswith("ENDATA"):
+                            break
+
+                        o = mosek_bas_re.match(line)
+                        if o is not None:
+                            if o.group(1) is not None:
+                                key = o.group(1)
+                                try:
+                                    skx[m.getvarnameindex(o.group(2))] = (
+                                        mosek.stakey.basis
+                                    )
+                                except:
+                                    pass
+                                try:
+                                    skc[m.getvarnameindex(o.group(3))] = (
+                                        mosek.stakey.low if key == "XL" else "XU"
+                                    )
+                                except:
+                                    pass
+                            else:
+                                key = o.group(4)
+                                name = o.group(5)
+                                stakey = (
+                                    mosek.stakey.low
+                                    if key == "LL"
+                                    else (
+                                        mosek.stakey.upr
+                                        if key == "UL"
+                                        else mosek.stakey.bas
+                                    )
+                                )
+
+                                try:
+                                    skx[m.getvarnameindex(name)] = stakey
+                                except:
+                                    try:
+                                        skc[m.getvarnameindex(name)] = stakey
+                                    except:
+                                        pass
+                m.putskc(mosek.soltype.bas, skc)
+                m.putskx(mosek.soltype.bas, skx)
             m.optimize()
 
             m.solutionsummary(mosek.streamtype.log)
 
             if basis_fn:
-                try:
-                    m.writedata(basis_fn)
-                except mosek.Error as err:
-                    logger.info("No model basis stored. Raised error: %s", err)
+                if m.solutiondef(mosek.soltype.bas):
+                    with open(basis_fn, "wt") as f:
+                        f.write(f"NAME {basis_fn}\n")
+
+                        skc = [
+                            (0 if sk != mosek.stakey.bas else 1, i, sk)
+                            for (i, sk) in enumerate(m.getskc(mosek.soltype.bas))
+                        ]
+                        skx = [
+                            (0 if sk == mosek.stakey.bas else 1, j, sk)
+                            for (j, sk) in enumerate(m.getskx(mosek.soltype.bas))
+                        ]
+                        skc.sort()
+                        skc.reverse()
+                        skx.sort()
+                        skx.reverse()
+                        numcon = m.getnumcon()
+                        while skx and skc and skx[-1][0] == 0 and skc[-1][0] == 0:
+                            (_, i, kc) = skc.pop()
+                            (_, j, kx) = skx.pop()
+
+                            namex = m.getvarname(j)
+                            namec = m.getconname(i)
+
+                            if kc in [mosek.stakey.low, mosek.stakey.fix]:
+                                f.write(f" XL {namex} {namec}\n")
+                            else:
+                                f.write(f" XU {namex} {namec}\n")
+                        while skc and skc[-1][0] == 0:
+                            (_, i, kc) = skc.pop()
+                            namec = m.getconname(i)
+                            if kc in [mosek.stakey.low, mosek.stakey.fix]:
+                                f.write(f" LL {namex}\n")
+                            else:
+                                f.write(f" UL {namex}\n")
+                        while skx:
+                            (_, j, kx) = skx.pop()
+                            namex = m.getvarname(j)
+                            if kx == mosek.stakey.bas:
+                                f.write(f" BS {namex}\n")
+                            elif kx in [mosek.stakey.low, mosek.stakey.fix]:
+                                f.write(f" LL {namex}\n")
+                            elif kx == mosek.stakey.upr:
+                                f.write(f" UL {namex}\n")
+                        f.write(f"ENDATA\n")
 
             soltype = None
             possible_soltypes = [

test/test_optimization.py

@@ -27,6 +27,13 @@
 if "highs" in available_solvers:
     params.append(("highs", "direct"))
 
+if "mosek" in available_solvers:
+    params.append(("mosek", "direct"))
+    params.append(("mosek", "lp"))
+elif "mosek_remote" in available_solvers:
+    params.append(("mosek_remote", "direct"))
+    params.append(("mosek_remote", "lp"))
+
 feasible_quadratic_solvers = quadratic_solvers
 # There seems to be a bug in scipopt with quadratic models on windows, see
 # https://github.com/PyPSA/linopy/actions/runs/7615240686/job/20739454099?pr=78
@@ -371,6 +378,10 @@ def test_solver_options(model, solver, io_api):
         "highs": {"time_limit": 1},
         "scip": {"limits/time": 1},
         "mosek": {"MSK_DPAR_OPTIMIZER_MAX_TIME": 1},
+        "mosek_remote": {
+            "MSK_DPAR_OPTIMIZER_MAX_TIME": 1,
+            "MSK_SPAR_REMOTE_OPTSERVER_HOST": "http://solve.mosek.com:30080",
+        },
         "mindopt": {"MaxTime": 1},
         "copt": {"TimeLimit": 1},
     }