diff --git a/linopy/io.py b/linopy/io.py index 823425c6..d7b9454c 100644 --- a/linopy/io.py +++ b/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_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) + + 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. diff --git a/linopy/model.py b/linopy/model.py index 06e9ed4a..2b604cc5 100644 --- a/linopy/model.py +++ b/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 diff --git a/linopy/solvers.py b/linopy/solvers.py index 99c2e5cf..e5aa00c0 100644 --- a/linopy/solvers.py +++ b/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,7 +870,9 @@ 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/ @@ -873,14 +887,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" ) 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 +905,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_mosekpy(m) + else: + m.readdata(problem_fn) for k, v in solver_options.items(): m.putparam(k, str(v)) @@ -899,17 +917,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 = [ diff --git a/test/test_optimization.py b/test/test_optimization.py index 98e157aa..18195982 100644 --- a/test/test_optimization.py +++ b/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 @@ -371,6 +379,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}, }