Coverage for picos/solvers/solver_mosek.py : 89.56%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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 (Optimizer)", "MOSEK via 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 (multidimensional) variable to a MOSEK scalar variable
135 (start) index."""
137 self._mosekBarVarIndex = dict()
138 """Maps suited PICOS symmetric variables to a MOSEK bar var index."""
140 self._mosekLinConOffset = dict()
141 """Maps a PICOS (multidimensional) linear constraint to a MOSEK scalar
142 constraint (start) index."""
144 self._mosekQuadConIndex = dict()
145 """Maps a PICOS quadratic constraint to a MOSEK constraint index."""
147 self._mosekCone = dict()
148 """
149 Maps a PICOS SOCC or RSOCC to a triple:
151 - The first entry is the index of a MOSEK conic constraint.
152 - The second entry is a list of MOSEK scalar (auxiliary) variable
153 indices that appear in the MOSEK conic constraint.
154 - The third entry is another list of same size containing auxiliary
155 constraints (or None). If an entry in the second list is None, then
156 the corresponding index in the first list is of a proper MOSEK scalar
157 variable instead of an auxiliary one, which can happen at most once
158 per variable and only if the structure of the PICOS constraint allows
159 it (that is, if the repspective entry of the cone member happens to be
160 a PICOS scalar variable as opposed to a composed affine expression).
161 """
163 self._mosekLMI = dict()
164 """Maps a PICOS LMI to a pair representing its MOSEK representation. The
165 first entry is the index of a MOSEK symmetric PSD "bar variable", the
166 second entry is the offset for a number of scalar linear equalities."""
168 self._mosekBarUnitCoefs = dict()
169 """Maps the row (column) count of a symmetric matrix to a list of MOSEK
170 symmetric coefficient matrices. More precisely, if X is a n×n symmetric
171 matrix, 0 ≤ k ≤ n*(n+1)/2, then <_mosekBarUnitCoefs[n][k], X> = h(X)[k]
172 where h refers to the lower-triangular, column-major half-vectorization
173 of X. These matrices are used to represent LMIs. Unlike values of
174 _mosekLMI, they are shared between LMIs of same size."""
176 def reset_problem(self):
177 """Implement :meth:`~.solver.Solver.reset_problem`."""
178 self.int = None
180 self._mosekVarOffset.clear()
181 self._mosekBarVarIndex.clear()
182 self._mosekLinConOffset.clear()
183 self._mosekQuadConIndex.clear()
184 self._mosekCone.clear()
185 self._mosekLMI.clear()
186 self._mosekBarUnitCoefs.clear()
188 @classmethod
189 def _get_environment(cls):
190 if not hasattr(cls, "mosekEnvironment"):
191 import mosek
192 cls.mosekEnvironment = mosek.Env()
194 return cls.mosekEnvironment
196 env = property(lambda self: self.__class__._get_environment())
197 """This references a MOSEK environment, which is shared among all
198 MOSEKSolver instances. (The MOSEK documentation states that "[a]ll tasks in
199 the program should share the same environment.")"""
201 @classmethod
202 def _get_major_version(cls):
203 if not hasattr(cls, "mosekVersion"):
204 import mosek
205 cls.mosekVersion = mosek.Env.getversion()
207 return cls.mosekVersion[0]
209 ver = property(lambda self: self.__class__._get_major_version())
210 """The major version of the available MOSEK library."""
212 @staticmethod
213 def _streamprinter(text):
214 sys.stdout.write(text)
215 sys.stdout.flush()
217 @staticmethod
218 def _low_tri_indices(rowCount):
219 """Yield enumerated lower triangular indices in column-major order."""
220 lti = 0
221 for col in range(rowCount):
222 for row in range(col, rowCount):
223 yield lti, row, col
224 lti += 1
226 def _index_function(self, picosVar, localIndex):
227 """Support :meth:`_affinexp_pic2msk`."""
228 if picosVar in self._mosekBarVarIndex:
229 assert picosVar not in self._mosekVarOffset
230 return None
231 else:
232 return self._mosekVarOffset[picosVar] + localIndex
234 def _affinexp_pic2msk(self, picosExpression):
235 """Convert a PICOS affine expression to a MOSEK one.
237 :yield:
238 A sequence of quintuples ``(J, V, K, W, c)``, each corresponding to
239 one scalar entry of a multidimensional PICOS affine expression,
240 where ``J`` are MOSEK scalar variable indices, ``V`` are associated
241 coefficients, ``K`` are MOSEK bar var (symmetric positive
242 semidefinite variable) indices, ``W`` are indices of associated
243 coefficient matrices and where ``c`` is a scalar constant.
244 """
245 # Handle the base case (no bar vars in the problem) quickly.
246 if not self._mosekBarVarIndex:
247 for J, V, c in picosExpression.sparse_rows(self._mosekVarOffset):
248 yield J, V, (), (), c
250 return
252 barVarIndices = tuple([] for _ in range(len(picosExpression)))
253 barVarCoefs = tuple([] for _ in range(len(picosExpression)))
255 for var, coef in picosExpression._linear_coefs.items():
256 if var in self._mosekBarVarIndex:
257 mosekBarVarIndex = self._mosekBarVarIndex[var]
259 assert isinstance(var, SymmetricVariable)
261 n = var.shape[0]
263 for localIndex in range(coef.size[0]):
264 localCoef = coef[localIndex, :]
266 if not localCoef:
267 continue
269 # Devectorize local coefficient.
270 # This works as c.T * svec(X) = <desvec(c), X> due to svec
271 # and its inverse desvec being isometric isomorphisms.
272 localCoefMatrix = var._vec.devectorize(localCoef.T)
274 assert localCoefMatrix.size == (n, n)
276 # Get lower triangular sparse rows for the coefficient.
277 # NOTE: This is much faster than using _low_low_tri_indices.
278 I, J, V = [], [], []
279 for j in range(n):
280 for i in range(j, n):
281 v = localCoefMatrix[i, j]
282 if v:
283 I.append(i)
284 J.append(j)
285 V.append(v)
287 # Create a new MOSEK symmetric coefficient matrix.
288 # NOTE: It does not seem possible to remove these matrices
289 # in MOSEK 9.2 so we might leak memory whenever they
290 # fall out of use (e.g. objective function updates).
291 mosekCoefMatrix = self.int.appendsparsesymmat(n, I, J, V)
293 # Append bar var index and coefficient matrix index.
294 barVarIndices[localIndex].append(mosekBarVarIndex)
295 barVarCoefs[localIndex].append(mosekCoefMatrix)
297 for localIndex, (J, V, c) in enumerate(picosExpression.sparse_rows(
298 varOffsetMap=None, indexFunction=self._index_function)):
299 # Remove indices and coefficients corresponding to bar vars.
300 JV = tuple(zip(*(jv for jv in zip(J, V) if jv[0] is not None)))
302 J, V = JV if JV else ((),)*2
303 K, W = barVarIndices[localIndex], barVarCoefs[localIndex]
305 yield J, V, K, W, c
307 def _scalar_affinexp_pic2msk(self, picosExpression):
308 assert len(picosExpression) == 1
309 return next(self._affinexp_pic2msk(picosExpression))
311 def _quadexp_pic2msk(self, picosExpression):
312 """Transform a quadratic expression from PICOS to MOSEK.
314 Tranforms the quadratic part of a PICOS quadratic expression to a
315 symmetric, sparse biliniar form of which only the lower triangular
316 entries are given, and that can be used with MOSEK's variable vector.
317 Note that MOSEK applies a built-in factor of 0.5 to all biliniar forms
318 while PICOS doesn't, so a factor of 2 is applied here to cancel it out.
319 """
320 assert isinstance(picosExpression, QuadraticExpression)
321 numVars = self.int.getnumvar()
323 # Mixing of conic and quadratic constraints is not supported by MOSEK;
324 # therefor we do not handle bar vars here.
325 assert not self._mosekBarVarIndex
327 # Make a sparse representation of the strict lower, diagonal and strict
328 # upper parts of the matrix.
329 IL, JL, VL = [], [], []
330 ID, JD, VD = [], [], []
331 IU, JU, VU = [], [], []
332 for (picosVar1, picosVar2), picosCoefs \
333 in picosExpression._sparse_quads.items():
334 for sparseIndex in range(len(picosCoefs)):
335 localVar1Index = picosCoefs.I[sparseIndex]
336 localVar2Index = picosCoefs.J[sparseIndex]
337 localCoefficient = picosCoefs.V[sparseIndex]
338 mskVar1Index = self._mosekVarOffset[picosVar1] + localVar1Index
339 mskVar2Index = self._mosekVarOffset[picosVar2] + localVar2Index
341 if mskVar2Index < mskVar1Index:
342 I, J, V = IL, JL, VL
343 elif mskVar1Index < mskVar2Index:
344 I, J, V = IU, JU, VU
345 else:
346 I, J, V = ID, JD, VD
348 I.append(mskVar1Index)
349 J.append(mskVar2Index)
350 V.append(localCoefficient)
352 # Compute the lower triangular part of the biliniar form.
353 L = cvxopt.spmatrix(VL, IL, JL, (numVars, numVars))
354 D = cvxopt.spmatrix(VD, ID, JD, (numVars, numVars))
355 U = cvxopt.spmatrix(VU, IU, JU, (numVars, numVars))
356 Q = 2*D + L + U.T
358 # Return it as a sparse triple for MOSEK to consume.
359 return list(Q.I), list(Q.J), list(Q.V)
361 def _import_variable(self, picosVar):
362 import mosek
364 numVars = self._mosekVarOffset[picosVar] = self.int.getnumvar()
365 dim = picosVar.dim
366 indices = range(numVars, numVars + dim)
367 self.int.appendvars(dim)
369 # Set the variable type.
370 if isinstance(picosVar, (IntegerVariable, BinaryVariable)):
371 self.int.putvartypelist(
372 indices, [mosek.variabletype.type_int]*dim)
374 # Import bounds, including the implied bounds of a BinaryVariable.
375 if isinstance(picosVar, BinaryVariable):
376 self.int.putvarboundlist(
377 indices, [mosek.boundkey.ra]*dim, [0]*dim, [1]*dim)
378 else:
379 boundKeys = [mosek.boundkey.fr]*dim
380 lowerBounds = [0.0]*dim
381 upperBounds = [0.0]*dim
383 lower, upper = picosVar.bound_dicts
385 for i, b in lower.items():
386 boundKeys[i] = mosek.boundkey.lo
387 lowerBounds[i] = b
389 for i, b in upper.items():
390 if boundKeys[i] == mosek.boundkey.fr:
391 boundKeys[i] = mosek.boundkey.up
392 else: # Also has a lower bound.
393 if lowerBounds[i] == b:
394 boundKeys[i] = mosek.boundkey.fx
395 else:
396 boundKeys[i] = mosek.boundkey.ra
398 upperBounds[i] = b
400 self.int.putvarboundlist(
401 indices, boundKeys, lowerBounds, upperBounds)
403 def _import_psd_variable(self, picosVar):
404 assert isinstance(picosVar, SymmetricVariable)
406 self._mosekBarVarIndex[picosVar] = self.int.getnumbarvar()
407 self.int.appendbarvars([picosVar.shape[0]])
409 def _import_linear_constraint(self, picosConstraint):
410 import mosek
412 numCons = self.int.getnumcon()
413 conLen = len(picosConstraint)
414 self.int.appendcons(conLen)
416 if picosConstraint.is_equality():
417 boundKey = mosek.boundkey.fx
418 elif picosConstraint.is_increasing():
419 boundKey = mosek.boundkey.up
420 else:
421 boundKey = mosek.boundkey.lo
423 rows = self._affinexp_pic2msk(picosConstraint.lhs - picosConstraint.rhs)
425 for localConIndex, (J, V, K, W, c) in enumerate(rows):
426 mosekConIndex = numCons + localConIndex
427 rhs = -c
429 self.int.putarow(mosekConIndex, J, V)
430 for k, w in zip(K, W):
431 self.int.putbaraij(mosekConIndex, k, [w], [1.0])
432 self.int.putconbound(mosekConIndex, boundKey, rhs, rhs)
434 self._mosekLinConOffset[picosConstraint] = numCons
436 def _import_quad_constraint(self, picosConstraint):
437 # Mixing of conic and quadratic constraints is not supported by MOSEK;
438 # therefor we do not handle bar vars here.
439 assert not self._mosekBarVarIndex
441 # Import the linear part first.
442 picosLinConPart = picosConstraint.le0.aff < 0
443 self._import_linear_constraint(picosLinConPart)
444 mosekConIndex = self._mosekLinConOffset.pop(picosLinConPart)
446 # Add the quadratic part.
447 self.int.putqconk(
448 mosekConIndex, *self._quadexp_pic2msk(picosConstraint.le0))
450 self._mosekQuadConIndex[picosConstraint] = mosekConIndex
452 def _var_was_used_in_cone(self, mosekVariableIndex, usedJustNow=[]):
453 if mosekVariableIndex in usedJustNow:
454 return True
455 for _, mosekVarIndices, _ in self._mosekCone.values():
456 if mosekVariableIndex in mosekVarIndices:
457 return True
458 return False
460 def _import_quad_conic_constraint(self, picosConstraint):
461 import mosek
463 isRotated = isinstance(picosConstraint, RSOCConstraint)
464 mosekVars, mosekCons = [], [None]*len(picosConstraint)
466 # Get an initial MOSEK representation of the cone member.
467 entries = []
468 if isRotated:
469 # MOSEK internally adds a factor of 2 to the upper bound while PICOS
470 # doesn't, so cancel it out by adding a factor of 1/sqrt(2) to both
471 # factors of the upper bound.
472 f = 1.0 / math.sqrt(2.0)
473 entries.append(self._scalar_affinexp_pic2msk(f*picosConstraint.ub1))
474 entries.append(self._scalar_affinexp_pic2msk(f*picosConstraint.ub2))
475 else:
476 entries.append(self._scalar_affinexp_pic2msk(picosConstraint.ub))
477 entries.extend(self._affinexp_pic2msk(picosConstraint.ne))
479 # Map cone member entries to existing MOSEK variables, if possible.
480 mosekVarsMissing = []
481 for scalarVarNum, (J, V, K, W, c) in enumerate(entries):
482 if len(J) == 1 and V[0] == 1.0 and not K and not W and not c \
483 and not self._var_was_used_in_cone(J[0], mosekVars):
484 mosekVars.append(J[0])
485 else:
486 mosekVars.append(None)
487 mosekVarsMissing.append(scalarVarNum)
489 # Create auxiliary variables and constraints.
490 numAux = len(mosekVarsMissing)
491 auxVarOffset = self.int.getnumvar()
492 auxConOffset = self.int.getnumcon()
493 self.int.appendvars(numAux)
494 self.int.appendcons(numAux)
496 # Mosek fixes (!) new variables at zero, so set them free.
497 self.int.putvarboundlist(range(auxVarOffset, auxVarOffset + numAux),
498 [mosek.boundkey.fr]*numAux, [0.0]*numAux, [0.0]*numAux)
500 # Constrain the auxiliary variables to be equal to the cone member
501 # entries for which no existing MOSEK variable could be used.
502 # TODO: Instead of always creating a constraint, fix variables via their
503 # bound whenever possible.
504 for auxNum, missingVarIndex in enumerate(mosekVarsMissing):
505 auxVarIndex = auxVarOffset + auxNum
506 auxConIndex = auxConOffset + auxNum
508 # Prepare the auxiliary constraint.
509 J, V, K, W, c = entries[missingVarIndex]
511 if self._debug():
512 self._debug(" Adding MOSEK auxiliary constraint: "
513 "{}.T * x{}{} = {}"
514 .format(V, J, " + f(barvars)" if K else "", -c))
516 # Add the auxiliary constraint.
517 self.int.putarow(auxConIndex, J, V)
518 self.int.putaij(auxConIndex, auxVarIndex, -1.0)
519 for k, w in zip(K, W):
520 self.int.putbaraij(auxConIndex, k, [w], [1.0])
521 self.int.putconbound(auxConIndex, mosek.boundkey.fx, -c, -c)
523 # Complete the mapping of cone member entries to MOSEK (auxiliary)
524 # variables (and auxiliary constraints).
525 mosekVars[missingVarIndex] = auxVarIndex
526 mosekCons[missingVarIndex] = auxConIndex
528 if self._debug():
529 self._debug(" Adding MOSEK conic constraint: {} in {}".format(
530 mosekVars, "Qr" if isRotated else "Q"))
532 # Add the conic constraint.
533 coneIndex = self.int.getnumcone()
534 mosekCone = mosek.conetype.rquad if isRotated else mosek.conetype.quad
535 self.int.appendcone(mosekCone, 0.0, mosekVars)
537 self._mosekCone[picosConstraint] = (coneIndex, mosekVars, mosekCons)
539 def _import_sdp_constraint(self, picosConstraint):
540 # NOTE: Trivial LMIs of the form X ≽ 0 are loaded as bar vars instead.
542 import mosek
544 n = picosConstraint.size[0]
545 dim = (n*(n + 1)) // 2
547 # MOSEK does not support general LMIs but so called "bar vars" which
548 # are variables in the symmetric positive semidefinite cone. We use them
549 # in combination with linear equalities to represent the LMI.
550 barVar = self.int.getnumbarvar()
551 mosekConOffset = self.int.getnumcon()
552 self.int.appendbarvars([n])
553 self.int.appendcons(dim)
555 # MOSEK uses a storage of symmetric coefficient matrices that are used
556 # as dot product coefficients to build scalar constraints involving both
557 # "bar vars" and normal scalar variables. We build a couple of these
558 # matrices to be able to select individual entries of our "bar vars".
559 # More precisely, if X is a n×n symmetric matrix and 0 ≤ k ≤ n*(n+1)/2,
560 # then <Units[n][k],X> = h(X)[k] where h refers to the lower-triangular,
561 # column-major half-vectorization of X.
562 if n in self._mosekBarUnitCoefs:
563 Units = self._mosekBarUnitCoefs[n]
564 else:
565 Units = self._mosekBarUnitCoefs[n] = [
566 self.int.appendsparsesymmat(
567 n, [row], [col], [1.0 if row == col else 0.5])
568 for col in range(n) for row in range(col, n)]
570 # We iterate over the lower triangular scalar sub-expressions of the
571 # expression that the PICOS constraint states to be PSD, and constrain
572 # them to be eqal to the MOSEK "bar var" at the same index.
573 psdRows = tuple(self._affinexp_pic2msk(picosConstraint.psd))
574 for lowTriIndex, row, col in self._low_tri_indices(n):
575 mosekConIndex = mosekConOffset + lowTriIndex
576 J, V, K, W, c = psdRows[row + n*col]
577 rhs = -c
579 # The lower-triangular entries in the PSD-constrained matrix …
580 self.int.putarow(mosekConIndex, J, V)
581 for k, w in zip(K, W):
582 self.int.putbaraij(mosekConIndex, k, [w], [1.0])
584 # … minus the corresponding bar var entries …
585 self.int.putbaraij(
586 mosekConIndex, barVar, [Units[lowTriIndex]], [-1.0])
588 # … should equal zero.
589 self.int.putconbound(mosekConIndex, mosek.boundkey.fx, rhs, rhs)
591 if self._debug():
592 self._debug(" Index {} ({}, {}): J = {}, V = {}, ..."
593 .format(lowTriIndex, row, col, J, V))
595 self._mosekLMI[picosConstraint] = (barVar, mosekConOffset)
597 def _import_constraint(self, picosConstraint):
598 if self._debug():
599 self._debug("Importing Constraint: {}".format(picosConstraint))
601 if isinstance(picosConstraint, AffineConstraint):
602 self._import_linear_constraint(picosConstraint)
603 elif isinstance(picosConstraint, ConvexQuadraticConstraint):
604 self._import_quad_constraint(picosConstraint)
605 elif isinstance(picosConstraint, SOCConstraint) \
606 or isinstance(picosConstraint, RSOCConstraint):
607 self._import_quad_conic_constraint(picosConstraint)
608 elif isinstance(picosConstraint, LMIConstraint):
609 self._import_sdp_constraint(picosConstraint)
610 else:
611 assert isinstance(picosConstraint, DummyConstraint), \
612 "Unexpected constraint type: {}".format(
613 picosConstraint.__class__.__name__)
615 def _reset_objective(self):
616 numVars = self.int.getnumvar()
617 numBarVars = self.int.getnumbarvar()
619 # Reset affine part.
620 self.int.putclist(range(numVars), [0.0]*numVars)
621 for k in range(numBarVars):
622 self.int.putbarcj(k, [], [])
623 self.int.putcfix(0.0)
625 # Reset quadratic part.
626 self.int.putqobj([], [], [])
628 def _import_affine_objective(self, picosObjective):
629 J, V, K, W, c = self._scalar_affinexp_pic2msk(picosObjective)
631 self.int.putclist(J, V)
632 for k, w in zip(K, W):
633 self.int.putbarcj(k, [w], [1.0])
634 self.int.putcfix(c)
636 def _import_quadratic_objective(self, picosObjective):
637 # Mixing of conic and quadratic constraints is not supported by MOSEK;
638 # therefor we do not handle bar vars here.
639 assert not self._mosekBarVarIndex
641 # Import the quadratic part.
642 self.int.putqobj(*self._quadexp_pic2msk(picosObjective))
644 # Import the affine part.
645 self._import_affine_objective(picosObjective.aff)
647 def _import_objective(self):
648 import mosek
650 picosSense, picosObjective = self.ext.no
652 # Import objective sense.
653 if picosSense == "min":
654 self.int.putobjsense(mosek.objsense.minimize)
655 else:
656 assert picosSense == "max"
657 self.int.putobjsense(mosek.objsense.maximize)
659 # Import objective function.
660 if isinstance(picosObjective, AffineExpression):
661 self._import_affine_objective(picosObjective)
662 else:
663 assert isinstance(picosObjective, QuadraticExpression)
664 self._import_quadratic_objective(picosObjective)
666 def _import_problem(self):
667 # Create a problem instance.
668 self.int = self.env.Task()
670 # Convert trivial LMIs and their symmetric variable into a PSD variable.
671 # This allows the constraint to be loaded more efficiently.
672 constraints, psdVars = [], set()
673 for constraint in self.ext.constraints.values():
674 if isinstance(constraint, LMIConstraint) and constraint.semidefVar:
675 psdVars.add(constraint.semidefVar)
676 else:
677 constraints.append(constraint)
679 # Import normal variables.
680 for variable in self.ext.variables.values():
681 if variable not in psdVars:
682 self._import_variable(variable)
684 # Import virtual PSD variables as MOSEK "bar variables".
685 for variable in psdVars:
686 self._import_psd_variable(variable)
688 # Import remaining constraints.
689 for constraint in constraints:
690 self._import_constraint(constraint)
692 # Set objective.
693 self._import_objective()
695 def _update_problem(self):
696 for oldConstraint in self._removed_constraints():
697 raise ProblemUpdateError("PICOS does not support removing "
698 "constraints from a MOSEK instance.")
700 for oldVariable in self._removed_variables():
701 raise ProblemUpdateError("PICOS does not support removing variables"
702 " from a MOSEK instance.")
704 for newVariable in self._new_variables():
705 self._import_variable(newVariable)
707 for newConstraint in self._new_constraints():
708 self._import_constraint(newConstraint)
710 if self._objective_has_changed():
711 self._reset_objective()
712 self._import_objective()
714 def _solve(self):
715 import mosek
717 # Determine whether a basic solution will be available.
718 # TODO: Give Problem an interface for checks like this.
719 isLP = isinstance(self.ext.no.function, AffineExpression) \
720 and all(isinstance(constraint, AffineConstraint)
721 for constraint in self.ext.constraints.values())
723 # Reset all solver options to default.
724 self.int.setdefaults()
725 self.int.putoptserverhost("")
727 # verbosity
728 if self._verbose():
729 self.int.set_Stream(mosek.streamtype.log, self._streamprinter)
730 self.int.putintparam(mosek.iparam.log, self.ext.verbosity())
732 # abs_prim_fsb_tol
733 if self.ext.options.abs_prim_fsb_tol is not None:
734 value = self.ext.options.abs_prim_fsb_tol
736 # Interior-point primal feasibility tolerances.
737 self.int.putdouparam(mosek.dparam.intpnt_tol_pfeas, value)
738 self.int.putdouparam(mosek.dparam.intpnt_co_tol_pfeas, value)
739 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_pfeas, value)
740 if self.ver <= 8:
741 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_pfeas, value)
743 # Simplex primal feasibility tolerance.
744 self.int.putdouparam(mosek.dparam.basis_tol_x, value)
746 # Mixed-integer (primal) feasibility tolerance.
747 self.int.putdouparam(mosek.dparam.mio_tol_feas, value)
749 # abs_dual_fsb_tol
750 if self.ext.options.abs_dual_fsb_tol is not None:
751 value = self.ext.options.abs_dual_fsb_tol
753 # Interior-point dual feasibility tolerances.
754 self.int.putdouparam(mosek.dparam.intpnt_tol_dfeas, value)
755 self.int.putdouparam(mosek.dparam.intpnt_co_tol_dfeas, value)
756 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_dfeas, value)
757 if self.ver <= 8:
758 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_dfeas, value)
760 # Simplex dual feasibility (optimality) tolerance.
761 self.int.putdouparam(mosek.dparam.basis_tol_s, value)
763 # rel_dual_fsb_tol
764 if self.ext.options.rel_dual_fsb_tol is not None:
765 # Simplex relative dual feasibility (optimality) tolerance.
766 self.int.putdouparam(mosek.dparam.basis_rel_tol_s,
767 self.ext.options.rel_dual_fsb_tol)
769 # rel_ipm_opt_tol
770 if self.ext.options.rel_ipm_opt_tol is not None:
771 value = self.ext.options.rel_ipm_opt_tol
773 # Interior-point primal feasibility tolerances.
774 self.int.putdouparam(mosek.dparam.intpnt_tol_rel_gap, value)
775 self.int.putdouparam(mosek.dparam.intpnt_co_tol_rel_gap, value)
776 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_rel_gap, value)
777 if self.ver <= 8:
778 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_rel_gap, value)
780 # abs_bnb_opt_tol
781 if self.ext.options.abs_bnb_opt_tol is not None:
782 self.int.putdouparam(mosek.dparam.mio_tol_abs_gap,
783 self.ext.options.abs_bnb_opt_tol)
785 # rel_bnb_opt_tol
786 if self.ext.options.rel_bnb_opt_tol is not None:
787 self.int.putdouparam(mosek.dparam.mio_tol_rel_gap,
788 self.ext.options.rel_bnb_opt_tol)
790 # integrality_tol
791 if self.ext.options.integrality_tol is not None:
792 self.int.putdouparam(mosek.dparam.mio_tol_abs_relax_int,
793 self.ext.options.integrality_tol)
795 # max_iterations
796 if self.ext.options.max_iterations is not None:
797 value = self.ext.options.max_iterations
798 self.int.putintparam(mosek.iparam.bi_max_iterations, value)
799 self.int.putintparam(mosek.iparam.intpnt_max_iterations, value)
800 self.int.putintparam(mosek.iparam.sim_max_iterations, value)
802 _lpm = {
803 "interior": (mosek.optimizertype.intpnt if isLP
804 else mosek.optimizertype.conic),
805 "psimplex": mosek.optimizertype.primal_simplex,
806 "dsimplex": mosek.optimizertype.dual_simplex}
808 # lp_node_method
809 if self.ext.options.lp_node_method is not None:
810 value = self.ext.options.lp_node_method
811 assert value in _lpm, "Unexpected lp_node_method value."
812 self.int.putintparam(mosek.iparam.mio_node_optimizer, _lpm[value])
814 # lp_root_method
815 if self.ext.options.lp_root_method is not None:
816 value = self.ext.options.lp_root_method
817 assert value in _lpm, "Unexpected lp_root_method value."
818 self.int.putintparam(mosek.iparam.mio_root_optimizer, _lpm[value])
820 # timelimit
821 if self.ext.options.timelimit is not None:
822 value = float(self.ext.options.timelimit)
823 self.int.putdouparam(mosek.dparam.optimizer_max_time, value)
824 self.int.putdouparam(mosek.dparam.mio_max_time, value)
826 # max_fsb_nodes
827 if self.ext.options.max_fsb_nodes is not None:
828 self.int.putintparam(mosek.iparam.mio_max_num_solutions,
829 self.ext.options.max_fsb_nodes)
831 # Handle MOSEK-specific options.
832 for key, value in self.ext.options.mosek_params.items():
833 try:
834 self.int.putparam(key.upper(), str(value))
835 except mosek.Error as error:
836 self._handle_bad_solver_specific_option(key, value, error)
838 # Handle 'mosek_server' option.
839 if self.ext.options.mosek_server:
840 self.int.putoptserverhost(self.ext.options.mosek_server)
842 # Handle unsupported options.
843 # TODO: Handle "hotstart" option (via mio_construct_sol).
844 self._handle_unsupported_option("hotstart", "treememory")
846 # Attempt to solve the problem.
847 with self._header(), self._stopwatch():
848 try:
849 self.int.optimize()
850 except mosek.Error as error:
851 if error.errno in (
852 mosek.rescode.err_con_q_not_psd,
853 mosek.rescode.err_con_q_not_nsd,
854 mosek.rescode.err_obj_q_not_psd,
855 mosek.rescode.err_obj_q_not_nsd,
856 mosek.rescode.err_toconic_constr_q_not_psd,
857 mosek.rescode.err_toconic_objective_not_psd):
858 self._handle_continuous_nonconvex_error(error)
859 else:
860 raise
862 # Set the solution to be retrieved.
863 if self.ext.is_continuous():
864 if isLP:
865 solType = mosek.soltype.bas
866 else:
867 solType = mosek.soltype.itr
868 else:
869 solType = mosek.soltype.itg
871 # Retrieve primals.
872 primals = {}
873 if self.ext.options.primals is not False:
874 values = [float("nan")]*self.int.getnumvar()
875 self.int.getxx(solType, values)
877 for picosVar in self.ext.variables.values():
878 if picosVar in self._mosekBarVarIndex:
879 barVarIndex = self._mosekBarVarIndex[picosVar]
881 assert isinstance(picosVar, SymmetricVariable)
883 n, d = picosVar.shape[0], picosVar.dim
884 lowTriPrimal = [float("nan")]*d
885 matrixPrimal = [float("nan")]*n**2
887 # Obtain a lower triangular vectorization from MOSEK.
888 self.int.getbarxj(solType, barVarIndex, lowTriPrimal)
890 # Convert to a full vectorization.
891 # TODO: Immediately convert to a symmetric vectorization.
892 for lti, row, col in self._low_tri_indices(n):
893 value = lowTriPrimal[lti]
894 matrixPrimal[n*row + col] = value
895 if row != col:
896 matrixPrimal[n*col + row] = value
898 # Load the full vectorization as a CVXOPT matrix.
899 matrixPrimal = cvxopt.matrix(matrixPrimal, (n, n))
901 # Re-vectorize using a symmetric vectorization.
902 primal = picosVar._vec.vectorize(matrixPrimal)
903 else:
904 mosekOffset = self._mosekVarOffset[picosVar]
905 primal = values[mosekOffset:mosekOffset + picosVar.dim]
907 if float("nan") in primal:
908 primals[picosVar] = None
909 else:
910 primals[picosVar] = primal
912 # Retrieve duals.
913 duals = {}
914 if self.ext.options.duals is not False and self.ext.is_continuous():
915 minimizing = self.int.getobjsense() == mosek.objsense.minimize
917 for constraint in self.ext.constraints.values():
918 if isinstance(constraint, DummyConstraint):
919 duals[constraint] = cvxopt.spmatrix(
920 [], [], [], constraint.size)
921 continue
923 length = len(constraint)
924 dual = [float("nan")]*length
926 if isinstance(constraint, AffineConstraint):
927 offset = self._mosekLinConOffset[constraint]
928 self.int.getyslice(solType, offset, offset + length, dual)
929 elif isinstance(constraint, ConvexQuadraticConstraint):
930 # TODO: Implement consistent QCQP dual retrieval for all
931 # solvers that return such duals.
932 dual = None
933 elif isinstance(constraint, SOCConstraint) \
934 or isinstance(constraint, RSOCConstraint):
935 mosekVars = self._mosekCone[constraint][1]
936 for localConeIndex in range(length):
937 x = [float("nan")]
938 offset = mosekVars[localConeIndex]
939 self.int.getsnxslice(solType, offset, offset + 1, x)
940 dual[localConeIndex] = x[0]
942 if isinstance(constraint, SOCConstraint):
943 dual = [-du for du in dual]
944 elif isinstance(constraint, RSOCConstraint):
945 dual[0] = -dual[0] / math.sqrt(2.0)
946 dual[1] = -dual[1] / math.sqrt(2.0)
947 dual[2:] = [-du for du in dual[2:]]
948 elif isinstance(constraint, LMIConstraint):
949 if constraint.semidefVar:
950 assert constraint not in self._mosekLMI
951 picosVar = constraint.semidefVar
952 assert picosVar in self._mosekBarVarIndex
953 barVar = self._mosekBarVarIndex[picosVar]
954 else:
955 assert not constraint.semidefVar
956 assert constraint in self._mosekLMI
957 barVar, _ = self._mosekLMI[constraint]
959 n = constraint.size[0]
960 d = n*(n + 1) // 2
961 lowerTriangularDual = [float("nan")]*d
962 self.int.getbarsj(solType, barVar, lowerTriangularDual)
963 for lti, row, col in self._low_tri_indices(n):
964 value = lowerTriangularDual[lti]
965 dual[n*row + col] = value
966 if row != col:
967 dual[n*col + row] = value
968 else:
969 assert False, "Constraint type not supported."
971 if dual is None:
972 pass
973 elif float("nan") in dual:
974 dual = None
975 else:
976 dual = cvxopt.matrix(dual, constraint.size)
978 if isinstance(constraint, AffineConstraint) \
979 or isinstance(constraint, LMIConstraint):
980 if not constraint.is_increasing():
981 dual = -dual
983 if minimizing:
984 dual = -dual
986 duals[constraint] = dual
988 # Retrieve objective value.
989 value = self.int.getprimalobj(solType)
991 # Retrieve solution status.
992 primalStatus = self._solution_status_pic2msk(
993 self.int.getsolsta(solType), primalSolution=True)
994 dualStatus = self._solution_status_pic2msk(
995 self.int.getsolsta(solType), primalSolution=False)
996 problemStatus = self._problem_status_pic2msk(
997 self.int.getprosta(solType), not self.ext.is_continuous())
999 return self._make_solution(
1000 value, primals, duals, primalStatus, dualStatus, problemStatus)
1002 def _solution_status_pic2msk(self, statusCode, primalSolution):
1003 from mosek import solsta as ss
1004 dualSolution = not primalSolution
1006 map = {
1007 ss.unknown: SS_UNKNOWN,
1008 ss.optimal: SS_OPTIMAL,
1009 ss.prim_feas:
1010 SS_FEASIBLE if primalSolution else SS_UNKNOWN,
1011 ss.dual_feas:
1012 SS_FEASIBLE if dualSolution else SS_UNKNOWN,
1013 ss.prim_and_dual_feas: SS_FEASIBLE,
1014 ss.prim_infeas_cer:
1015 SS_INFEASIBLE if primalSolution else SS_UNKNOWN,
1016 ss.dual_infeas_cer:
1017 SS_INFEASIBLE if dualSolution else SS_UNKNOWN,
1018 ss.prim_illposed_cer: SS_UNKNOWN,
1019 ss.dual_illposed_cer: SS_UNKNOWN,
1020 ss.integer_optimal: SS_OPTIMAL,
1021 }
1023 if self.ver < 9:
1024 map.update({
1025 ss.near_optimal: SS_FEASIBLE,
1026 ss.near_prim_feas: SS_UNKNOWN,
1027 ss.near_dual_feas: SS_UNKNOWN,
1028 ss.near_prim_and_dual_feas: SS_UNKNOWN,
1029 ss.near_prim_infeas_cer: SS_UNKNOWN,
1030 ss.near_dual_infeas_cer: SS_UNKNOWN,
1031 ss.near_integer_optimal: SS_FEASIBLE
1032 })
1034 try:
1035 return map[statusCode]
1036 except KeyError:
1037 self._warn(
1038 "The MOSEK solution status code {} is not known to PICOS."
1039 .format(statusCode))
1040 return SS_UNKNOWN
1042 def _problem_status_pic2msk(self, statusCode, integerProblem):
1043 from mosek import prosta as ps
1045 map = {
1046 ps.unknown: PS_UNKNOWN,
1047 ps.prim_and_dual_feas: PS_FEASIBLE,
1048 ps.prim_feas:
1049 PS_FEASIBLE if integerProblem else PS_UNKNOWN,
1050 ps.dual_feas: PS_UNKNOWN,
1051 ps.prim_infeas: PS_INFEASIBLE,
1052 ps.dual_infeas: PS_UNBOUNDED, # TODO: UNB_OR_INF
1053 ps.prim_and_dual_infeas: PS_INFEASIBLE,
1054 ps.ill_posed: PS_ILLPOSED,
1055 ps.prim_infeas_or_unbounded: PS_UNKNOWN # TODO: UNB_OR_INF
1056 }
1058 if self.ver < 9:
1059 map.update({
1060 ps.near_prim_and_dual_feas: PS_UNKNOWN,
1061 ps.near_prim_feas: PS_UNKNOWN,
1062 ps.near_dual_feas: PS_UNKNOWN
1063 })
1065 try:
1066 return map[statusCode]
1067 except KeyError:
1068 self._warn("The MOSEK problem status code {} is not known to PICOS."
1069 .format(statusCode))
1070 return PS_UNKNOWN
1073# --------------------------------------
1074__all__ = api_end(_API_START, globals())