Coverage for picos/solvers/solver_mosek.py: 88.61%
553 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2018-2019 Maximilian Stahlberg
3#
4# This file is part of PICOS.
5#
6# PICOS is free software: you can redistribute it and/or modify it under the
7# terms of the GNU General Public License as published by the Free Software
8# Foundation, either version 3 of the License, or (at your option) any later
9# version.
10#
11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
17# ------------------------------------------------------------------------------
19"""Implementation of :class:`MOSEKSolver`."""
21import math
22import sys
24import cvxopt
26from ..apidoc import api_end, api_start
27from ..constraints import (AffineConstraint, ConvexQuadraticConstraint,
28 DummyConstraint, LMIConstraint, RSOCConstraint,
29 SOCConstraint)
30from ..expressions import (AffineExpression, BinaryVariable, IntegerVariable,
31 QuadraticExpression, SymmetricVariable)
32from ..modeling.footprint import Specification
33from ..modeling.solution import (PS_FEASIBLE, PS_ILLPOSED, PS_INFEASIBLE,
34 PS_UNBOUNDED, PS_UNKNOWN, SS_FEASIBLE,
35 SS_INFEASIBLE, SS_OPTIMAL, SS_UNKNOWN)
36from .solver import ProblemUpdateError, Solver
38_API_START = api_start(globals())
39# -------------------------------
42class MOSEKSolver(Solver):
43 """Interface to the MOSEK solver via its low level Optimizer API.
45 Supports both MOSEK 8 and 9.
47 The low level API is tedious to interface, but is currently much faster than
48 the high level Fusion API, which would be the prefered interface otherwise.
49 """
51 SUPPORTED = Specification(
52 objectives=[
53 AffineExpression,
54 QuadraticExpression],
55 constraints=[
56 DummyConstraint,
57 AffineConstraint,
58 SOCConstraint,
59 RSOCConstraint,
60 ConvexQuadraticConstraint,
61 LMIConstraint])
63 @classmethod
64 def supports(cls, footprint, explain=False):
65 """Implement :meth:`~.solver.Solver.supports`."""
66 result = Solver.supports(footprint, explain)
67 if not result or (explain and not result[0]):
68 return result
70 # No nonconvex quadratic objectives.
71 if footprint.nonconvex_quadratic_objective:
72 if explain:
73 return False, "Problems with nonconvex quadratic objectives."
74 else:
75 return False
77 conic = any(("con", constraint) in footprint
78 for constraint in (SOCConstraint, RSOCConstraint, LMIConstraint))
80 quadratic = footprint.objective.clstype is QuadraticExpression \
81 or ("con", ConvexQuadraticConstraint) in footprint
83 # No mixing of quadratic and conic problems.
84 if conic and quadratic:
85 if explain:
86 return (False,
87 "Problems that mix conic and quadratic constraints.")
88 else:
89 return False
91 # No integer SDPs.
92 if footprint.integer and ("con", LMIConstraint) in footprint:
93 if explain:
94 return False, "Integer Semidefinite Programs."
95 else:
96 return False
98 if footprint not in cls.SUPPORTED:
99 if explain:
100 return False, cls.SUPPORTED.mismatch_reason(footprint)
101 else:
102 return False
104 return (True, None) if explain else True
106 @classmethod
107 def default_penalty(cls):
108 """Implement :meth:`~.solver.Solver.default_penalty`."""
109 return 0.0 # Commercial solver.
111 @classmethod
112 def test_availability(cls):
113 """Implement :meth:`~.solver.Solver.test_availability`."""
114 cls.check_import("mosek")
116 @classmethod
117 def names(cls):
118 """Implement :meth:`~.solver.Solver.names`."""
119 return "mosek", "MOSEK", "MOSEK", "Optimizer API"
121 @classmethod
122 def is_free(cls):
123 """Implement :meth:`~.solver.Solver.is_free`."""
124 return False
126 def __init__(self, problem):
127 """Initialize a MOSEK (Optimizer) solver interface.
129 :param ~picos.Problem problem: The problem to be solved.
130 """
131 super(MOSEKSolver, self).__init__(problem)
133 self._mosekVarOffset = dict()
134 """Maps a PICOS variable to a MOSEK scalar variable offset."""
136 self._mosekBarVarIndex = dict()
137 """Maps suited PICOS symmetric variables to a MOSEK bar var index."""
139 self._mosekLinConOffset = dict()
140 """Maps a PICOS linear constraint to a MOSEK constraint offset."""
142 self._mosekQuadConIndex = dict()
143 """Maps a PICOS quadratic constraint to a MOSEK constraint index."""
145 self._mosekCone = dict()
146 """
147 Maps a PICOS SOCC or RSOCC to a triple:
149 - The first entry is the index of a MOSEK conic constraint.
150 - The second entry is a list of MOSEK scalar (auxiliary) variable
151 indices that appear in the MOSEK conic constraint.
152 - The third entry is another list of same size containing auxiliary
153 constraints (or None). If an entry in the second list is None, then
154 the corresponding index in the first list is of a proper MOSEK scalar
155 variable instead of an auxiliary one, which can happen at most once
156 per variable and only if the structure of the PICOS constraint allows
157 it (that is, if the repspective entry of the cone member happens to be
158 a PICOS scalar variable as opposed to a composed affine expression).
159 """
161 self._mosekLMI = dict()
162 """Maps a PICOS LMI to a pair representing its MOSEK representation. The
163 first entry is the index of a MOSEK symmetric PSD "bar variable", the
164 second entry is the offset for a number of scalar linear equalities."""
166 self._mosekBarUnitCoefs = dict()
167 """Maps the row (column) count of a symmetric matrix to a list of MOSEK
168 symmetric coefficient matrices. More precisely, if X is a n×n symmetric
169 matrix, 0 ≤ k ≤ n*(n+1)/2, then <_mosekBarUnitCoefs[n][k], X> = h(X)[k]
170 where h refers to the lower-triangular, column-major half-vectorization
171 of X. These matrices are used to represent LMIs. Unlike values of
172 _mosekLMI, they are shared between LMIs of same size."""
174 def reset_problem(self):
175 """Implement :meth:`~.solver.Solver.reset_problem`."""
176 self.int = None
178 self._mosekVarOffset.clear()
179 self._mosekBarVarIndex.clear()
180 self._mosekLinConOffset.clear()
181 self._mosekQuadConIndex.clear()
182 self._mosekCone.clear()
183 self._mosekLMI.clear()
184 self._mosekBarUnitCoefs.clear()
186 @classmethod
187 def _get_environment(cls):
188 if not hasattr(cls, "mosekEnvironment"):
189 import mosek
190 cls.mosekEnvironment = mosek.Env()
192 return cls.mosekEnvironment
194 env = property(lambda self: self.__class__._get_environment())
195 """This references a MOSEK environment, which is shared among all
196 MOSEKSolver instances. (The MOSEK documentation states that "[a]ll tasks in
197 the program should share the same environment.")"""
199 @classmethod
200 def _get_major_version(cls):
201 if not hasattr(cls, "mosekVersion"):
202 import mosek
203 cls.mosekVersion = mosek.Env.getversion()
205 return cls.mosekVersion[0]
207 ver = property(lambda self: self.__class__._get_major_version())
208 """The major version of the available MOSEK library."""
210 @staticmethod
211 def _streamprinter(text):
212 sys.stdout.write(text)
213 sys.stdout.flush()
215 @staticmethod
216 def _low_tri_indices(rowCount):
217 """Yield enumerated lower triangular indices in column-major order."""
218 lti = 0
219 for col in range(rowCount):
220 for row in range(col, rowCount):
221 yield lti, row, col
222 lti += 1
224 def _affinexp_pic2msk(self, picosExpression):
225 """Convert a PICOS affine expression to a MOSEK one.
227 :yield:
228 A sequence of quintuples ``(J, V, K, W, c)``, each corresponding to
229 one scalar entry of a multidimensional PICOS affine expression,
230 where ``J`` are MOSEK scalar variable indices, ``V`` are associated
231 coefficients, ``K`` are MOSEK bar var (symmetric positive
232 semidefinite variable) indices, ``W`` are indices of associated
233 coefficient matrices and where ``c`` is a scalar constant.
234 """
235 # Handle the base case (no bar vars in the problem) quickly.
236 if not self._mosekBarVarIndex:
237 for J, V, c in picosExpression.sparse_rows(self._mosekVarOffset):
238 yield J, V, (), (), c
240 return
242 barVarIndices = tuple([] for _ in range(len(picosExpression)))
243 barVarCoefs = tuple([] for _ in range(len(picosExpression)))
245 for var, coef in picosExpression._linear_coefs.items():
246 if var in self._mosekBarVarIndex:
247 mosekBarVarIndex = self._mosekBarVarIndex[var]
249 assert isinstance(var, SymmetricVariable)
251 n = var.shape[0]
253 for localIndex in range(coef.size[0]):
254 localCoef = coef[localIndex, :]
256 if not localCoef:
257 continue
259 # Devectorize local coefficient.
260 # This works as c.T * svec(X) = <desvec(c), X> due to svec
261 # and its inverse desvec being isometric isomorphisms.
262 localCoefMatrix = var._vec.devectorize(localCoef.T)
264 assert localCoefMatrix.size == (n, n)
266 # Get lower triangular sparse rows for the coefficient.
267 # NOTE: This is much faster than using _low_low_tri_indices.
268 I, J, V = [], [], []
269 for j in range(n):
270 for i in range(j, n):
271 v = localCoefMatrix[i, j]
272 if v:
273 I.append(i)
274 J.append(j)
275 V.append(v)
277 # Create a new MOSEK symmetric coefficient matrix.
278 # NOTE: It does not seem possible to remove these matrices
279 # in MOSEK 9.2 so we might leak memory whenever they
280 # fall out of use (e.g. objective function updates).
281 mosekCoefMatrix = self.int.appendsparsesymmat(n, I, J, V)
283 # Append bar var index and coefficient matrix index.
284 barVarIndices[localIndex].append(mosekBarVarIndex)
285 barVarCoefs[localIndex].append(mosekCoefMatrix)
287 rows = enumerate(picosExpression.sparse_rows(self._mosekVarOffset))
288 for localIndex, (J, V, c) in rows:
289 K, W = barVarIndices[localIndex], barVarCoefs[localIndex]
291 yield J, V, K, W, c
293 def _scalar_affinexp_pic2msk(self, picosExpression):
294 assert len(picosExpression) == 1
295 return next(self._affinexp_pic2msk(picosExpression))
297 def _quadexp_pic2msk(self, picosExpression):
298 """Transform a quadratic expression from PICOS to MOSEK.
300 Tranforms the quadratic part of a PICOS quadratic expression to a
301 symmetric, sparse biliniar form of which only the lower triangular
302 entries are given, and that can be used with MOSEK's variable vector.
303 Note that MOSEK applies a built-in factor of 0.5 to all biliniar forms
304 while PICOS doesn't, so a factor of 2 is applied here to cancel it out.
305 """
306 assert isinstance(picosExpression, QuadraticExpression)
307 numVars = self.int.getnumvar()
309 # Mixing of conic and quadratic constraints is not supported by MOSEK;
310 # therefor we do not handle bar vars here.
311 assert not self._mosekBarVarIndex
313 # Make a sparse representation of the strict lower, diagonal and strict
314 # upper parts of the matrix.
315 IL, JL, VL = [], [], []
316 ID, JD, VD = [], [], []
317 IU, JU, VU = [], [], []
318 for (picosVar1, picosVar2), picosCoefs \
319 in picosExpression._sparse_quads.items():
320 for sparseIndex in range(len(picosCoefs)):
321 localVar1Index = picosCoefs.I[sparseIndex]
322 localVar2Index = picosCoefs.J[sparseIndex]
323 localCoefficient = picosCoefs.V[sparseIndex]
324 mskVar1Index = self._mosekVarOffset[picosVar1] + localVar1Index
325 mskVar2Index = self._mosekVarOffset[picosVar2] + localVar2Index
327 if mskVar2Index < mskVar1Index:
328 I, J, V = IL, JL, VL
329 elif mskVar1Index < mskVar2Index:
330 I, J, V = IU, JU, VU
331 else:
332 I, J, V = ID, JD, VD
334 I.append(mskVar1Index)
335 J.append(mskVar2Index)
336 V.append(localCoefficient)
338 # Compute the lower triangular part of the biliniar form.
339 L = cvxopt.spmatrix(VL, IL, JL, (numVars, numVars))
340 D = cvxopt.spmatrix(VD, ID, JD, (numVars, numVars))
341 U = cvxopt.spmatrix(VU, IU, JU, (numVars, numVars))
342 Q = 2*D + L + U.T
344 # Return it as a sparse triple for MOSEK to consume.
345 return list(Q.I), list(Q.J), list(Q.V)
347 def _import_variable(self, picosVar):
348 import mosek
350 numVars = self._mosekVarOffset[picosVar] = self.int.getnumvar()
351 dim = picosVar.dim
352 indices = range(numVars, numVars + dim)
353 self.int.appendvars(dim)
355 # Set the variable type.
356 if isinstance(picosVar, (IntegerVariable, BinaryVariable)):
357 self.int.putvartypelist(
358 indices, [mosek.variabletype.type_int]*dim)
360 # Import bounds, including the implied bounds of a BinaryVariable.
361 if isinstance(picosVar, BinaryVariable):
362 self.int.putvarboundlist(
363 indices, [mosek.boundkey.ra]*dim, [0]*dim, [1]*dim)
364 else:
365 boundKeys = [mosek.boundkey.fr]*dim
366 lowerBounds = [0.0]*dim
367 upperBounds = [0.0]*dim
369 lower, upper = picosVar.bound_dicts
371 for i, b in lower.items():
372 boundKeys[i] = mosek.boundkey.lo
373 lowerBounds[i] = b
375 for i, b in upper.items():
376 if boundKeys[i] == mosek.boundkey.fr:
377 boundKeys[i] = mosek.boundkey.up
378 else: # Also has a lower bound.
379 if lowerBounds[i] == b:
380 boundKeys[i] = mosek.boundkey.fx
381 else:
382 boundKeys[i] = mosek.boundkey.ra
384 upperBounds[i] = b
386 self.int.putvarboundlist(
387 indices, boundKeys, lowerBounds, upperBounds)
389 def _import_psd_variable(self, picosVar):
390 assert isinstance(picosVar, SymmetricVariable)
392 self._mosekBarVarIndex[picosVar] = self.int.getnumbarvar()
393 self.int.appendbarvars([picosVar.shape[0]])
395 def _import_linear_constraint(self, picosConstraint):
396 import mosek
398 numCons = self.int.getnumcon()
399 conLen = len(picosConstraint)
400 self.int.appendcons(conLen)
402 if picosConstraint.is_equality():
403 boundKey = mosek.boundkey.fx
404 elif picosConstraint.is_increasing():
405 boundKey = mosek.boundkey.up
406 else:
407 boundKey = mosek.boundkey.lo
409 rows = self._affinexp_pic2msk(picosConstraint.lhs - picosConstraint.rhs)
411 for localConIndex, (J, V, K, W, c) in enumerate(rows):
412 mosekConIndex = numCons + localConIndex
413 rhs = -c
415 self.int.putarow(mosekConIndex, J, V)
416 for k, w in zip(K, W):
417 self.int.putbaraij(mosekConIndex, k, [w], [1.0])
418 self.int.putconbound(mosekConIndex, boundKey, rhs, rhs)
420 self._mosekLinConOffset[picosConstraint] = numCons
422 def _import_quad_constraint(self, picosConstraint):
423 # Mixing of conic and quadratic constraints is not supported by MOSEK;
424 # therefor we do not handle bar vars here.
425 assert not self._mosekBarVarIndex
427 # Import the linear part first.
428 picosLinConPart = picosConstraint.le0.aff < 0
429 self._import_linear_constraint(picosLinConPart)
430 mosekConIndex = self._mosekLinConOffset.pop(picosLinConPart)
432 # Add the quadratic part.
433 self.int.putqconk(
434 mosekConIndex, *self._quadexp_pic2msk(picosConstraint.le0))
436 self._mosekQuadConIndex[picosConstraint] = mosekConIndex
438 def _var_was_used_in_cone(self, mosekVariableIndex, usedJustNow=[]):
439 if mosekVariableIndex in usedJustNow:
440 return True
441 for _, mosekVarIndices, _ in self._mosekCone.values():
442 if mosekVariableIndex in mosekVarIndices:
443 return True
444 return False
446 def _import_quad_conic_constraint(self, picosConstraint):
447 import mosek
449 isRotated = isinstance(picosConstraint, RSOCConstraint)
450 mosekVars, mosekCons = [], [None]*len(picosConstraint)
452 # Get an initial MOSEK representation of the cone member.
453 entries = []
454 if isRotated:
455 # MOSEK internally adds a factor of 2 to the upper bound while PICOS
456 # doesn't, so cancel it out by adding a factor of 1/sqrt(2) to both
457 # factors of the upper bound.
458 f = 1.0 / math.sqrt(2.0)
459 entries.append(self._scalar_affinexp_pic2msk(f*picosConstraint.ub1))
460 entries.append(self._scalar_affinexp_pic2msk(f*picosConstraint.ub2))
461 else:
462 entries.append(self._scalar_affinexp_pic2msk(picosConstraint.ub))
463 entries.extend(self._affinexp_pic2msk(picosConstraint.ne))
465 # Map cone member entries to existing MOSEK variables, if possible.
466 mosekVarsMissing = []
467 for scalarVarNum, (J, V, K, W, c) in enumerate(entries):
468 if len(J) == 1 and V[0] == 1.0 and not K and not W and not c \
469 and not self._var_was_used_in_cone(J[0], mosekVars):
470 mosekVars.append(J[0])
471 else:
472 mosekVars.append(None)
473 mosekVarsMissing.append(scalarVarNum)
475 # Create auxiliary variables and constraints.
476 numAux = len(mosekVarsMissing)
477 auxVarOffset = self.int.getnumvar()
478 auxConOffset = self.int.getnumcon()
479 self.int.appendvars(numAux)
480 self.int.appendcons(numAux)
482 # Mosek fixes (!) new variables at zero, so set them free.
483 self.int.putvarboundlist(range(auxVarOffset, auxVarOffset + numAux),
484 [mosek.boundkey.fr]*numAux, [0.0]*numAux, [0.0]*numAux)
486 # Constrain the auxiliary variables to be equal to the cone member
487 # entries for which no existing MOSEK variable could be used.
488 # TODO: Instead of always creating a constraint, fix variables via their
489 # bound whenever possible.
490 for auxNum, missingVarIndex in enumerate(mosekVarsMissing):
491 auxVarIndex = auxVarOffset + auxNum
492 auxConIndex = auxConOffset + auxNum
494 # Prepare the auxiliary constraint.
495 J, V, K, W, c = entries[missingVarIndex]
497 if self._debug():
498 self._debug(" Adding MOSEK auxiliary constraint: "
499 "{}.T * x{}{} = {}"
500 .format(V, J, " + f(barvars)" if K else "", -c))
502 # Add the auxiliary constraint.
503 self.int.putarow(auxConIndex, J, V)
504 self.int.putaij(auxConIndex, auxVarIndex, -1.0)
505 for k, w in zip(K, W):
506 self.int.putbaraij(auxConIndex, k, [w], [1.0])
507 self.int.putconbound(auxConIndex, mosek.boundkey.fx, -c, -c)
509 # Complete the mapping of cone member entries to MOSEK (auxiliary)
510 # variables (and auxiliary constraints).
511 mosekVars[missingVarIndex] = auxVarIndex
512 mosekCons[missingVarIndex] = auxConIndex
514 if self._debug():
515 self._debug(" Adding MOSEK conic constraint: {} in {}".format(
516 mosekVars, "Qr" if isRotated else "Q"))
518 # Add the conic constraint.
519 coneIndex = self.int.getnumcone()
520 mosekCone = mosek.conetype.rquad if isRotated else mosek.conetype.quad
521 self.int.appendcone(mosekCone, 0.0, mosekVars)
523 self._mosekCone[picosConstraint] = (coneIndex, mosekVars, mosekCons)
525 def _import_sdp_constraint(self, picosConstraint):
526 # NOTE: Trivial LMIs of the form X ≽ 0 are loaded as bar vars instead.
528 import mosek
530 n = picosConstraint.size[0]
531 dim = (n*(n + 1)) // 2
533 # MOSEK does not support general LMIs but so called "bar vars" which
534 # are variables in the symmetric positive semidefinite cone. We use them
535 # in combination with linear equalities to represent the LMI.
536 barVar = self.int.getnumbarvar()
537 mosekConOffset = self.int.getnumcon()
538 self.int.appendbarvars([n])
539 self.int.appendcons(dim)
541 # MOSEK uses a storage of symmetric coefficient matrices that are used
542 # as dot product coefficients to build scalar constraints involving both
543 # "bar vars" and normal scalar variables. We build a couple of these
544 # matrices to be able to select individual entries of our "bar vars".
545 # More precisely, if X is a n×n symmetric matrix and 0 ≤ k ≤ n*(n+1)/2,
546 # then <Units[n][k],X> = h(X)[k] where h refers to the lower-triangular,
547 # column-major half-vectorization of X.
548 if n in self._mosekBarUnitCoefs:
549 Units = self._mosekBarUnitCoefs[n]
550 else:
551 Units = self._mosekBarUnitCoefs[n] = [
552 self.int.appendsparsesymmat(
553 n, [row], [col], [1.0 if row == col else 0.5])
554 for col in range(n) for row in range(col, n)]
556 # We iterate over the lower triangular scalar sub-expressions of the
557 # expression that the PICOS constraint states to be PSD, and constrain
558 # them to be eqal to the MOSEK "bar var" at the same index.
559 psdRows = tuple(self._affinexp_pic2msk(picosConstraint.psd))
560 for lowTriIndex, row, col in self._low_tri_indices(n):
561 mosekConIndex = mosekConOffset + lowTriIndex
562 J, V, K, W, c = psdRows[row + n*col]
563 rhs = -c
565 # The lower-triangular entries in the PSD-constrained matrix …
566 self.int.putarow(mosekConIndex, J, V)
567 for k, w in zip(K, W):
568 self.int.putbaraij(mosekConIndex, k, [w], [1.0])
570 # … minus the corresponding bar var entries …
571 self.int.putbaraij(
572 mosekConIndex, barVar, [Units[lowTriIndex]], [-1.0])
574 # … should equal zero.
575 self.int.putconbound(mosekConIndex, mosek.boundkey.fx, rhs, rhs)
577 if self._debug():
578 self._debug(" Index {} ({}, {}): J = {}, V = {}, ..."
579 .format(lowTriIndex, row, col, J, V))
581 self._mosekLMI[picosConstraint] = (barVar, mosekConOffset)
583 def _import_constraint(self, picosConstraint):
584 if self._debug():
585 self._debug("Importing Constraint: {}".format(picosConstraint))
587 if isinstance(picosConstraint, AffineConstraint):
588 self._import_linear_constraint(picosConstraint)
589 elif isinstance(picosConstraint, ConvexQuadraticConstraint):
590 self._import_quad_constraint(picosConstraint)
591 elif isinstance(picosConstraint, SOCConstraint) \
592 or isinstance(picosConstraint, RSOCConstraint):
593 self._import_quad_conic_constraint(picosConstraint)
594 elif isinstance(picosConstraint, LMIConstraint):
595 self._import_sdp_constraint(picosConstraint)
596 else:
597 assert isinstance(picosConstraint, DummyConstraint), \
598 "Unexpected constraint type: {}".format(
599 picosConstraint.__class__.__name__)
601 def _reset_objective(self):
602 numVars = self.int.getnumvar()
603 numBarVars = self.int.getnumbarvar()
605 # Reset affine part.
606 self.int.putclist(range(numVars), [0.0]*numVars)
607 for k in range(numBarVars):
608 self.int.putbarcj(k, [], [])
609 self.int.putcfix(0.0)
611 # Reset quadratic part.
612 self.int.putqobj([], [], [])
614 def _import_affine_objective(self, picosObjective):
615 J, V, K, W, c = self._scalar_affinexp_pic2msk(picosObjective)
617 self.int.putclist(J, V)
618 for k, w in zip(K, W):
619 self.int.putbarcj(k, [w], [1.0])
620 self.int.putcfix(c)
622 def _import_quadratic_objective(self, picosObjective):
623 # Mixing of conic and quadratic constraints is not supported by MOSEK;
624 # therefor we do not handle bar vars here.
625 assert not self._mosekBarVarIndex
627 # Import the quadratic part.
628 self.int.putqobj(*self._quadexp_pic2msk(picosObjective))
630 # Import the affine part.
631 self._import_affine_objective(picosObjective.aff)
633 def _import_objective(self):
634 import mosek
636 picosSense, picosObjective = self.ext.no
638 # Import objective sense.
639 if picosSense == "min":
640 self.int.putobjsense(mosek.objsense.minimize)
641 else:
642 assert picosSense == "max"
643 self.int.putobjsense(mosek.objsense.maximize)
645 # Import objective function.
646 if isinstance(picosObjective, AffineExpression):
647 self._import_affine_objective(picosObjective)
648 else:
649 assert isinstance(picosObjective, QuadraticExpression)
650 self._import_quadratic_objective(picosObjective)
652 def _import_problem(self):
653 # Create a problem instance.
654 self.int = self.env.Task()
656 # Convert trivial LMIs and their symmetric variable into a PSD variable.
657 # This allows the constraint to be loaded more efficiently.
658 constraints, psdVars = [], set()
659 for constraint in self.ext.constraints.values():
660 if isinstance(constraint, LMIConstraint) and constraint.semidefVar:
661 psdVars.add(constraint.semidefVar)
662 else:
663 constraints.append(constraint)
665 # Import normal variables.
666 for variable in self.ext.variables.values():
667 if variable not in psdVars:
668 self._import_variable(variable)
670 # Import virtual PSD variables as MOSEK "bar variables".
671 for variable in psdVars:
672 self._import_psd_variable(variable)
674 # Import remaining constraints.
675 for constraint in constraints:
676 self._import_constraint(constraint)
678 # Set objective.
679 self._import_objective()
681 def _update_problem(self):
682 for oldConstraint in self._removed_constraints():
683 raise ProblemUpdateError("PICOS does not support removing "
684 "constraints from a MOSEK instance.")
686 for oldVariable in self._removed_variables():
687 raise ProblemUpdateError("PICOS does not support removing variables"
688 " from a MOSEK instance.")
690 for newVariable in self._new_variables():
691 self._import_variable(newVariable)
693 for newConstraint in self._new_constraints():
694 self._import_constraint(newConstraint)
696 if self._objective_has_changed():
697 self._reset_objective()
698 self._import_objective()
700 def _solve(self):
701 import mosek
703 # Determine whether an LP is being solved.
704 # TODO: Give Problem an interface for checks like this.
705 _affine = (AffineConstraint, DummyConstraint)
706 _is_lp = isinstance(self.ext.no.function, AffineExpression) \
707 and all(isinstance(constraint, _affine)
708 for constraint in self.ext.constraints.values())
710 # Reset all solver options to default.
711 self.int.setdefaults()
712 self.int.putoptserverhost("")
714 # verbosity
715 if self._verbose():
716 self.int.set_Stream(mosek.streamtype.log, self._streamprinter)
717 self.int.putintparam(mosek.iparam.log, self.ext.verbosity())
719 # abs_prim_fsb_tol
720 if self.ext.options.abs_prim_fsb_tol is not None:
721 value = self.ext.options.abs_prim_fsb_tol
723 # Interior-point primal feasibility tolerances.
724 self.int.putdouparam(mosek.dparam.intpnt_tol_pfeas, value)
725 self.int.putdouparam(mosek.dparam.intpnt_co_tol_pfeas, value)
726 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_pfeas, value)
727 if self.ver <= 8:
728 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_pfeas, value)
730 # Simplex primal feasibility tolerance.
731 self.int.putdouparam(mosek.dparam.basis_tol_x, value)
733 # Mixed-integer (primal) feasibility tolerance.
734 self.int.putdouparam(mosek.dparam.mio_tol_feas, value)
736 # abs_dual_fsb_tol
737 if self.ext.options.abs_dual_fsb_tol is not None:
738 value = self.ext.options.abs_dual_fsb_tol
740 # Interior-point dual feasibility tolerances.
741 self.int.putdouparam(mosek.dparam.intpnt_tol_dfeas, value)
742 self.int.putdouparam(mosek.dparam.intpnt_co_tol_dfeas, value)
743 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_dfeas, value)
744 if self.ver <= 8:
745 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_dfeas, value)
747 # Simplex dual feasibility (optimality) tolerance.
748 self.int.putdouparam(mosek.dparam.basis_tol_s, value)
750 # rel_dual_fsb_tol
751 if self.ext.options.rel_dual_fsb_tol is not None:
752 # Simplex relative dual feasibility (optimality) tolerance.
753 self.int.putdouparam(mosek.dparam.basis_rel_tol_s,
754 self.ext.options.rel_dual_fsb_tol)
756 # rel_ipm_opt_tol
757 if self.ext.options.rel_ipm_opt_tol is not None:
758 value = self.ext.options.rel_ipm_opt_tol
760 # Interior-point primal feasibility tolerances.
761 self.int.putdouparam(mosek.dparam.intpnt_tol_rel_gap, value)
762 self.int.putdouparam(mosek.dparam.intpnt_co_tol_rel_gap, value)
763 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_rel_gap, value)
764 if self.ver <= 8:
765 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_rel_gap, value)
767 # abs_bnb_opt_tol
768 if self.ext.options.abs_bnb_opt_tol is not None:
769 self.int.putdouparam(mosek.dparam.mio_tol_abs_gap,
770 self.ext.options.abs_bnb_opt_tol)
772 # rel_bnb_opt_tol
773 if self.ext.options.rel_bnb_opt_tol is not None:
774 self.int.putdouparam(mosek.dparam.mio_tol_rel_gap,
775 self.ext.options.rel_bnb_opt_tol)
777 # integrality_tol
778 if self.ext.options.integrality_tol is not None:
779 self.int.putdouparam(mosek.dparam.mio_tol_abs_relax_int,
780 self.ext.options.integrality_tol)
782 # max_iterations
783 if self.ext.options.max_iterations is not None:
784 value = self.ext.options.max_iterations
785 self.int.putintparam(mosek.iparam.bi_max_iterations, value)
786 self.int.putintparam(mosek.iparam.intpnt_max_iterations, value)
787 self.int.putintparam(mosek.iparam.sim_max_iterations, value)
789 # Prepare lp_node_method and lp_root_method.
790 _lpm = {
791 "interior": (mosek.optimizertype.intpnt if _is_lp
792 else mosek.optimizertype.conic),
793 "psimplex": mosek.optimizertype.primal_simplex,
794 "dsimplex": mosek.optimizertype.dual_simplex}
796 # lp_node_method
797 if self.ext.options.lp_node_method is not None:
798 value = self.ext.options.lp_node_method
799 assert value in _lpm, "Unexpected lp_node_method value."
800 self.int.putintparam(mosek.iparam.mio_node_optimizer, _lpm[value])
802 # lp_root_method
803 if self.ext.options.lp_root_method is not None:
804 value = self.ext.options.lp_root_method
805 assert value in _lpm, "Unexpected lp_root_method value."
806 self.int.putintparam(mosek.iparam.mio_root_optimizer, _lpm[value])
808 # timelimit
809 if self.ext.options.timelimit is not None:
810 value = float(self.ext.options.timelimit)
811 self.int.putdouparam(mosek.dparam.optimizer_max_time, value)
812 self.int.putdouparam(mosek.dparam.mio_max_time, value)
814 # max_fsb_nodes
815 if self.ext.options.max_fsb_nodes is not None:
816 self.int.putintparam(mosek.iparam.mio_max_num_solutions,
817 self.ext.options.max_fsb_nodes)
819 # Handle MOSEK-specific options.
820 for key, value in self.ext.options.mosek_params.items():
821 try:
822 self.int.putparam(key.upper(), str(value))
823 except mosek.Error as error:
824 self._handle_bad_solver_specific_option(key, value, error)
826 # Handle 'mosek_basic_sol' option.
827 if self.ext.options.mosek_basic_sol and _is_lp:
828 _intpnt_basis = mosek.basindtype.always
829 else:
830 _intpnt_basis = mosek.basindtype.never
831 self.int.putintparam(mosek.iparam.intpnt_basis, _intpnt_basis)
833 # Handle 'mosek_server' option.
834 if self.ext.options.mosek_server:
835 self.int.putoptserverhost(self.ext.options.mosek_server)
837 # Handle unsupported options.
838 # TODO: Handle "hotstart" option (via mio_construct_sol).
839 self._handle_unsupported_option("hotstart", "treememory")
841 # Attempt to solve the problem.
842 with self._header(), self._stopwatch():
843 try:
844 self.int.optimize()
845 except mosek.Error as error:
846 if error.errno in (
847 mosek.rescode.err_con_q_not_psd,
848 mosek.rescode.err_con_q_not_nsd,
849 mosek.rescode.err_obj_q_not_psd,
850 mosek.rescode.err_obj_q_not_nsd,
851 mosek.rescode.err_toconic_constr_q_not_psd,
852 mosek.rescode.err_toconic_objective_not_psd):
853 self._handle_continuous_nonconvex_error(error)
854 else:
855 raise
857 # Set the solution to be retrieved.
858 if self.ext.is_continuous():
859 if _intpnt_basis == mosek.basindtype.always:
860 solType = mosek.soltype.bas
861 else:
862 solType = mosek.soltype.itr
863 else:
864 solType = mosek.soltype.itg
866 # Retrieve primals.
867 primals = {}
868 if self.ext.options.primals is not False:
869 values = [float("nan")]*self.int.getnumvar()
870 self.int.getxx(solType, values)
872 for picosVar in self.ext.variables.values():
873 if picosVar in self._mosekBarVarIndex:
874 barVarIndex = self._mosekBarVarIndex[picosVar]
876 assert isinstance(picosVar, SymmetricVariable)
878 n, d = picosVar.shape[0], picosVar.dim
879 lowTriPrimal = [float("nan")]*d
880 matrixPrimal = [float("nan")]*n**2
882 # Obtain a lower triangular vectorization from MOSEK.
883 self.int.getbarxj(solType, barVarIndex, lowTriPrimal)
885 # Convert to a full vectorization.
886 # TODO: Immediately convert to a symmetric vectorization.
887 for lti, row, col in self._low_tri_indices(n):
888 value = lowTriPrimal[lti]
889 matrixPrimal[n*row + col] = value
890 if row != col:
891 matrixPrimal[n*col + row] = value
893 # Load the full vectorization as a CVXOPT matrix.
894 matrixPrimal = cvxopt.matrix(matrixPrimal, (n, n))
896 # Re-vectorize using a symmetric vectorization.
897 primal = picosVar._vec.vectorize(matrixPrimal)
898 else:
899 mosekOffset = self._mosekVarOffset[picosVar]
900 primal = values[mosekOffset:mosekOffset + picosVar.dim]
902 if float("nan") in primal:
903 primals[picosVar] = None
904 else:
905 primals[picosVar] = primal
907 # Retrieve duals.
908 duals = {}
909 if self.ext.options.duals is not False and self.ext.is_continuous():
910 minimizing = self.int.getobjsense() == mosek.objsense.minimize
912 for constraint in self.ext.constraints.values():
913 if isinstance(constraint, DummyConstraint):
914 duals[constraint] = cvxopt.spmatrix(
915 [], [], [], constraint.size)
916 continue
918 length = len(constraint)
919 dual = [float("nan")]*length
921 if isinstance(constraint, AffineConstraint):
922 offset = self._mosekLinConOffset[constraint]
923 self.int.getyslice(solType, offset, offset + length, dual)
924 elif isinstance(constraint, ConvexQuadraticConstraint):
925 # TODO: Implement consistent QCQP dual retrieval for all
926 # solvers that return such duals.
927 dual = None
928 elif isinstance(constraint, SOCConstraint) \
929 or isinstance(constraint, RSOCConstraint):
930 mosekVars = self._mosekCone[constraint][1]
931 for localConeIndex in range(length):
932 x = [float("nan")]
933 offset = mosekVars[localConeIndex]
934 self.int.getsnxslice(solType, offset, offset + 1, x)
935 dual[localConeIndex] = x[0]
937 if isinstance(constraint, SOCConstraint):
938 dual = [-du for du in dual]
939 elif isinstance(constraint, RSOCConstraint):
940 dual[0] = -dual[0] / math.sqrt(2.0)
941 dual[1] = -dual[1] / math.sqrt(2.0)
942 dual[2:] = [-du for du in dual[2:]]
943 elif isinstance(constraint, LMIConstraint):
944 if constraint.semidefVar:
945 assert constraint not in self._mosekLMI
946 picosVar = constraint.semidefVar
947 assert picosVar in self._mosekBarVarIndex
948 barVar = self._mosekBarVarIndex[picosVar]
949 else:
950 assert not constraint.semidefVar
951 assert constraint in self._mosekLMI
952 barVar, _ = self._mosekLMI[constraint]
954 n = constraint.size[0]
955 d = n*(n + 1) // 2
956 lowerTriangularDual = [float("nan")]*d
957 self.int.getbarsj(solType, barVar, lowerTriangularDual)
958 for lti, row, col in self._low_tri_indices(n):
959 value = lowerTriangularDual[lti]
960 dual[n*row + col] = value
961 if row != col:
962 dual[n*col + row] = value
963 else:
964 assert False, "Constraint type not supported."
966 if dual is None:
967 pass
968 elif float("nan") in dual:
969 dual = None
970 else:
971 dual = cvxopt.matrix(dual, constraint.size)
973 if isinstance(constraint, AffineConstraint) \
974 or isinstance(constraint, LMIConstraint):
975 if not constraint.is_increasing():
976 dual = -dual
978 if minimizing:
979 dual = -dual
981 duals[constraint] = dual
983 # Retrieve objective value.
984 value = self.int.getprimalobj(solType)
986 # Retrieve solution status.
987 primalStatus = self._solution_status_pic2msk(
988 self.int.getsolsta(solType), primalSolution=True)
989 dualStatus = self._solution_status_pic2msk(
990 self.int.getsolsta(solType), primalSolution=False)
991 problemStatus = self._problem_status_pic2msk(
992 self.int.getprosta(solType), not self.ext.is_continuous())
994 return self._make_solution(
995 value, primals, duals, primalStatus, dualStatus, problemStatus)
997 def _solution_status_pic2msk(self, statusCode, primalSolution):
998 from mosek import solsta as ss
999 dualSolution = not primalSolution
1001 map = {
1002 ss.unknown: SS_UNKNOWN,
1003 ss.optimal: SS_OPTIMAL,
1004 ss.prim_feas:
1005 SS_FEASIBLE if primalSolution else SS_UNKNOWN,
1006 ss.dual_feas:
1007 SS_FEASIBLE if dualSolution else SS_UNKNOWN,
1008 ss.prim_and_dual_feas: SS_FEASIBLE,
1009 ss.prim_infeas_cer:
1010 SS_INFEASIBLE if primalSolution else SS_UNKNOWN,
1011 ss.dual_infeas_cer:
1012 SS_INFEASIBLE if dualSolution else SS_UNKNOWN,
1013 ss.prim_illposed_cer: SS_UNKNOWN,
1014 ss.dual_illposed_cer: SS_UNKNOWN,
1015 ss.integer_optimal: SS_OPTIMAL,
1016 }
1018 if self.ver < 9:
1019 map.update({
1020 ss.near_optimal: SS_FEASIBLE,
1021 ss.near_prim_feas: SS_UNKNOWN,
1022 ss.near_dual_feas: SS_UNKNOWN,
1023 ss.near_prim_and_dual_feas: SS_UNKNOWN,
1024 ss.near_prim_infeas_cer: SS_UNKNOWN,
1025 ss.near_dual_infeas_cer: SS_UNKNOWN,
1026 ss.near_integer_optimal: SS_FEASIBLE
1027 })
1029 try:
1030 return map[statusCode]
1031 except KeyError:
1032 self._warn(
1033 "The MOSEK solution status code {} is not known to PICOS."
1034 .format(statusCode))
1035 return SS_UNKNOWN
1037 def _problem_status_pic2msk(self, statusCode, integerProblem):
1038 from mosek import prosta as ps
1040 map = {
1041 ps.unknown: PS_UNKNOWN,
1042 ps.prim_and_dual_feas: PS_FEASIBLE,
1043 ps.prim_feas:
1044 PS_FEASIBLE if integerProblem else PS_UNKNOWN,
1045 ps.dual_feas: PS_UNKNOWN,
1046 ps.prim_infeas: PS_INFEASIBLE,
1047 ps.dual_infeas: PS_UNBOUNDED, # TODO: UNB_OR_INF
1048 ps.prim_and_dual_infeas: PS_INFEASIBLE,
1049 ps.ill_posed: PS_ILLPOSED,
1050 ps.prim_infeas_or_unbounded: PS_UNKNOWN # TODO: UNB_OR_INF
1051 }
1053 if self.ver < 9:
1054 map.update({
1055 ps.near_prim_and_dual_feas: PS_UNKNOWN,
1056 ps.near_prim_feas: PS_UNKNOWN,
1057 ps.near_dual_feas: PS_UNKNOWN
1058 })
1060 try:
1061 return map[statusCode]
1062 except KeyError:
1063 self._warn("The MOSEK problem status code {} is not known to PICOS."
1064 .format(statusCode))
1065 return PS_UNKNOWN
1068# --------------------------------------
1069__all__ = api_end(_API_START, globals())