Merge pull request #244 from PyPSA/mosek-direct

Add quadratic objective support and 'direct' mode for MOSEK (superset of #226)

doc/release_notes.rst

@@ -4,6 +4,11 @@ Release Notes
 Upcoming Version
 ----------------
 
+
+**New Features**
+
+* A direct interface to the `Mosek` solver was added. With this change, a new conversion function `model.to_mosek` was added to convert a linopy model to a `mosek` model. The `solve` function now supports the `mosek` solver with `io_api="direct"`.
+
 * It is now possible to create LinearExpression from a `pandas.DataFrame`, `pandas.Series`, a `numpy.array` or constant scalar values, e.g. `linopy.LinearExpression(df)`. This will create a LinearExpression with constants only and the coordinates of the DataFrame, Series or array as dimensions.
 
 **Bugfixes**

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,104 @@ def to_file(m, fn, integer_label="general"):
     return fn
 
 
+def to_mosek(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)
+
+    labels = M.vlabels.astype(str).astype(object)
+    task.generatevarnames(
+        np.arange(0, len(labels)), "x%0", [len(labels)], None, [0], list(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 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_mosek,
+    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_mosek = to_mosek
+
     to_highspy = to_highspy
 
     to_block_files = to_block_files

linopy/solvers.py

@@ -77,6 +77,7 @@
             m.checkinall()
 
         available_solvers.append("mosek")
+
 with contextlib.suppress(ImportError):
     import mindoptpy
 
@@ -94,7 +95,15 @@
 
 
 io_structure = dict(
-    lp_file={"gurobi", "xpress", "cbc", "glpk", "cplex", "mosek", "mindopt"},
+    lp_file={
+        "gurobi",
+        "xpress",
+        "cbc",
+        "glpk",
+        "cplex",
+        "mosek",
+        "mindopt",
+    },
     blocks={"pips"},
 )
 
@@ -845,6 +854,9 @@ def get_solver_solution() -> Solution:
     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(
     model,
     io_api=None,
@@ -858,12 +870,20 @@ def run_mosek(
     **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/
 
     For more information on solver options, see
     https://docs.mosek.com/latest/pythonapi/parameters.html#doc-all-parameter-list
+
+
+    For remote optimization of smaller problems, which do not require a license,
+    set the following solver_options:
+    {"MSK_SPAR_REMOTE_OPTSERVER_HOST": "http://solve.mosek.com:30080"}
+
     """
     CONDITION_MAP = {
         "solsta.unknown": "unknown",
@@ -873,14 +893,15 @@ def run_mosek(
         "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"
+            "Keyword argument `io_api` has to be one of `lp`, `mps`, `direct` 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)
@@ -890,7 +911,10 @@ def run_mosek(
             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_mosek(m)
+            else:
+                m.readdata(problem_fn)
 
             for k, v in solver_options.items():
                 m.putparam(k, str(v))
@@ -899,17 +923,107 @@ def run_mosek(
                 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,14 @@
 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