Coverage for picos/solvers/solver_mosek.py : 12.55%
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)
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._mosekLinConOffset = dict()
138 """Maps a PICOS (multidimensional) linear constraint to a MOSEK scalar
139 constraint (start) index."""
141 self._mosekQuadConIndex = dict()
142 """Maps a PICOS quadratic constraint to a MOSEK constraint index."""
144 self._mosekCone = dict()
145 """
146 Maps a PICOS SOCC or RSOCC to a triple:
148 - The first entry is the index of a MOSEK conic constraint.
149 - The second entry is a list of MOSEK scalar (auxiliary) variable
150 indices that appear in the MOSEK conic constraint.
151 - The third entry is another list of same size containing auxiliary
152 constraints (or None). If an entry in the second list is None, then
153 the corresponding index in the first list is of a proper MOSEK scalar
154 variable instead of an auxiliary one, which can happen at most once
155 per variable and only if the structure of the PICOS constraint allows
156 it (that is, if the repspective entry of the cone member happens to be
157 a PICOS scalar variable as opposed to a composed affine expression).
158 """
160 self._mosekLMI = dict()
161 """Maps a PICOS LMI to a pair representing its MOSEK representation. The
162 first entry is the index of a MOSEK symmetric PSD "bar variable", the
163 second entry is the offset for a number of scalar linear equalities."""
165 self._mosekBarUnitCoefs = dict()
166 """Maps the row (column) count of a symmetric matrix to a list of MOSEK
167 symmetric coefficient matrices. More precisely, if X is a n×n symmetric
168 matrix, 0 ≤ k ≤ n*(n+1)/2, then <_mosekBarUnitCoefs[n][k], X> = h(X)[k]
169 where h refers to the lower-triangular, column-major half-vectorization
170 of X. These matrices are used to represent LMIs. Unlike values of
171 _mosekLMI, they are shared between LMIs of same size."""
173 def reset_problem(self):
174 """Implement :meth:`~.solver.Solver.reset_problem`."""
175 self.int = None
177 self._mosekVarOffset.clear()
178 self._mosekLinConOffset.clear()
179 self._mosekQuadConIndex.clear()
180 self._mosekCone.clear()
181 self._mosekLMI.clear()
182 self._mosekBarUnitCoefs.clear()
184 @classmethod
185 def _get_environment(cls):
186 if not hasattr(cls, "mosekEnvironment"):
187 import mosek
188 cls.mosekEnvironment = mosek.Env()
190 return cls.mosekEnvironment
192 env = property(lambda self: self.__class__._get_environment())
193 """This references a MOSEK environment, which is shared among all
194 MOSEKSolver instances. (The MOSEK documentation states that "[a]ll tasks in
195 the program should share the same environment.")"""
197 @classmethod
198 def _get_major_version(cls):
199 if not hasattr(cls, "mosekVersion"):
200 import mosek
201 cls.mosekVersion = mosek.Env.getversion()
203 return cls.mosekVersion[0]
205 ver = property(lambda self: self.__class__._get_major_version())
206 """The major version of the available MOSEK library."""
208 @staticmethod
209 def _streamprinter(text):
210 sys.stdout.write(text)
211 sys.stdout.flush()
213 @staticmethod
214 def _low_tri_indices(rowCount):
215 """Yield lower triangular (row, col) indices in column-major order."""
216 for col in range(rowCount):
217 for row in range(col, rowCount):
218 yield (row, col)
220 def _scalar_affinexp_pic2msk(self, picosExpression):
221 assert len(picosExpression) == 1
222 return picosExpression.sparse_rows(self._mosekVarOffset)[0]
224 def _affinexp_pic2msk(self, picosExpression):
225 return picosExpression.sparse_rows(self._mosekVarOffset)
227 def _quadexp_pic2msk(self, picosExpression):
228 """Transform a quadratic expression from PICOS to MOSEK.
230 Tranforms the quadratic part of a PICOS quadratic expression to a
231 symmetric, sparse biliniar form of which only the lower triangular
232 entries are given, and that can be used with MOSEK's variable vector.
233 Note that MOSEK applies a built-in factor of 0.5 to all biliniar forms
234 while PICOS doesn't, so a factor of 2 is applied here to cancel it out.
235 """
236 assert isinstance(picosExpression, QuadraticExpression)
237 numVars = self.int.getnumvar()
239 # Make a sparse representation of the strict lower, diagonal and strict
240 # upper parts of the matrix.
241 IL, JL, VL = [], [], []
242 ID, JD, VD = [], [], []
243 IU, JU, VU = [], [], []
244 for (picosVar1, picosVar2), picosCoefs \
245 in picosExpression._sparse_quads.items():
246 for sparseIndex in range(len(picosCoefs)):
247 localVar1Index = picosCoefs.I[sparseIndex]
248 localVar2Index = picosCoefs.J[sparseIndex]
249 localCoefficient = picosCoefs.V[sparseIndex]
250 mskVar1Index = self._mosekVarOffset[picosVar1] + localVar1Index
251 mskVar2Index = self._mosekVarOffset[picosVar2] + localVar2Index
253 if mskVar2Index < mskVar1Index:
254 I, J, V = IL, JL, VL
255 elif mskVar1Index < mskVar2Index:
256 I, J, V = IU, JU, VU
257 else:
258 I, J, V = ID, JD, VD
260 I.append(mskVar1Index)
261 J.append(mskVar2Index)
262 V.append(localCoefficient)
264 # Compute the lower triangular part of the biliniar form.
265 L = cvxopt.spmatrix(VL, IL, JL, (numVars, numVars))
266 D = cvxopt.spmatrix(VD, ID, JD, (numVars, numVars))
267 U = cvxopt.spmatrix(VU, IU, JU, (numVars, numVars))
268 Q = 2*D + L + U.T
270 # Return it as a sparse triple for MOSEK to consume.
271 return list(Q.I), list(Q.J), list(Q.V)
273 def _import_variable(self, picosVar):
274 import mosek
276 numVars = self._mosekVarOffset[picosVar] = self.int.getnumvar()
277 dim = picosVar.dim
278 indices = range(numVars, numVars + dim)
279 self.int.appendvars(dim)
281 # Set the variable type.
282 if isinstance(picosVar, (IntegerVariable, BinaryVariable)):
283 self.int.putvartypelist(
284 indices, [mosek.variabletype.type_int]*dim)
286 # Import bounds, including the implied bounds of a BinaryVariable.
287 if isinstance(picosVar, BinaryVariable):
288 self.int.putvarboundlist(
289 indices, [mosek.boundkey.ra]*dim, [0]*dim, [1]*dim)
290 else:
291 boundKeys = [mosek.boundkey.fr]*dim
292 lowerBounds = [0.0]*dim
293 upperBounds = [0.0]*dim
295 lower, upper = picosVar.bound_dicts
297 for i, b in lower.items():
298 boundKeys[i] = mosek.boundkey.lo
299 lowerBounds[i] = b
301 for i, b in upper.items():
302 if boundKeys[i] == mosek.boundkey.fr:
303 boundKeys[i] = mosek.boundkey.up
304 else: # Also has a lower bound.
305 if lowerBounds[i] == b:
306 boundKeys[i] = mosek.boundkey.fx
307 else:
308 boundKeys[i] = mosek.boundkey.ra
310 upperBounds[i] = b
312 self.int.putvarboundlist(
313 indices, boundKeys, lowerBounds, upperBounds)
315 def _import_linear_constraint(self, picosConstraint):
316 import mosek
318 numCons = self.int.getnumcon()
319 conLen = len(picosConstraint)
320 self.int.appendcons(conLen)
322 rows = picosConstraint.sparse_Ab_rows(self._mosekVarOffset)
324 if picosConstraint.is_equality():
325 boundKey = mosek.boundkey.fx
326 elif picosConstraint.is_increasing():
327 boundKey = mosek.boundkey.up
328 else:
329 boundKey = mosek.boundkey.lo
331 for localConIndex, (mosekVarIndices, coefs, offset) in enumerate(rows):
332 mosekConIndex = numCons + localConIndex
334 self.int.putarow(mosekConIndex, mosekVarIndices, coefs)
335 self.int.putconbound(mosekConIndex, boundKey, offset, offset)
337 self._mosekLinConOffset[picosConstraint] = numCons
339 def _import_quad_constraint(self, picosConstraint):
340 # Import the linear part first.
341 picosLinConPart = picosConstraint.le0.aff < 0
342 self._import_linear_constraint(picosLinConPart)
343 mosekConIndex = self._mosekLinConOffset.pop(picosLinConPart)
345 # Add the quadratic part.
346 self.int.putqconk(
347 mosekConIndex, *self._quadexp_pic2msk(picosConstraint.le0))
349 self._mosekQuadConIndex[picosConstraint] = mosekConIndex
351 def _var_was_used_in_cone(self, mosekVariableIndex, usedJustNow=[]):
352 if mosekVariableIndex in usedJustNow:
353 return True
354 for _, mosekVarIndices, _ in self._mosekCone.values():
355 if mosekVariableIndex in mosekVarIndices:
356 return True
357 return False
359 def _import_quad_conic_constraint(self, picosConstraint):
360 import mosek
362 isRotated = isinstance(picosConstraint, RSOCConstraint)
363 mosekVars, mosekCons = [], [None]*len(picosConstraint)
365 # Get an initial MOSEK representation of the cone member.
366 entries = []
367 if isRotated:
368 # MOSEK internally adds a factor of 2 to the upper bound while PICOS
369 # doesn't, so cancel it out by adding a factor of 1/sqrt(2) to both
370 # factors of the upper bound.
371 f = 1.0 / math.sqrt(2.0)
372 entries.append(self._scalar_affinexp_pic2msk(f*picosConstraint.ub1))
373 entries.append(self._scalar_affinexp_pic2msk(f*picosConstraint.ub2))
374 else:
375 entries.append(self._scalar_affinexp_pic2msk(picosConstraint.ub))
376 entries += self._affinexp_pic2msk(picosConstraint.ne)
378 # Map cone member entries to existing MOSEK variables, if possible.
379 mosekVarsMissing = []
380 for scalarVarNum, (indices, values, offset) in enumerate(entries):
381 if len(values) == 1 and values[0] == 1.0 and offset == 0.0 \
382 and not self._var_was_used_in_cone(indices[0], mosekVars):
383 mosekVars.append(indices[0])
384 else:
385 mosekVars.append(None)
386 mosekVarsMissing.append(scalarVarNum)
388 # Create auxiliary variables and constraints.
389 numAux = len(mosekVarsMissing)
390 auxVarOffset = self.int.getnumvar()
391 auxConOffset = self.int.getnumcon()
392 self.int.appendvars(numAux)
393 self.int.appendcons(numAux)
395 # Mosek fixes (!) new variables at zero, so set them free.
396 self.int.putvarboundlist(range(auxVarOffset, auxVarOffset + numAux),
397 [mosek.boundkey.fr]*numAux, [0.0]*numAux, [0.0]*numAux)
399 # Constrain the auxiliary variables to be equal to the cone member
400 # entries for which no existing MOSEK variable could be used.
401 for auxNum, missingVarIndex in enumerate(mosekVarsMissing):
402 auxVarIndex = auxVarOffset + auxNum
403 auxConIndex = auxConOffset + auxNum
405 # Prepare the auxiliary constraint.
406 indices, values, offset = entries[missingVarIndex]
408 # TODO: Instead of always creating a constraint, fix variables via
409 # their bound whenever possible.
410 # if len(indices) is 0:
411 # if self._debug():
412 # self._debug(" Fixing MOSEK auxiliary variable: "
413 # "x[{}] = {}".format(auxVarIndex, offset))
414 #
415 # self.int.putvarbound(
416 # auxVarIndex, mosek.boundkey.fx, offset, offset)
417 # else:
418 indices.append(auxVarIndex)
419 values.append(-1.0)
421 if self._debug():
422 self._debug(" Adding MOSEK auxiliary constraint: "
423 "{}.T * x{} = {}".format(values, indices, -offset))
425 # Add the auxiliary constraint.
426 self.int.putarow(auxConIndex, indices, values)
427 self.int.putconbound(
428 auxConIndex, mosek.boundkey.fx, -offset, -offset)
430 # Complete the mapping of cone member entries to MOSEK (auxiliary)
431 # variables (and auxiliary constraints).
432 mosekVars[missingVarIndex] = auxVarIndex
433 mosekCons[missingVarIndex] = auxConIndex
435 if self._debug():
436 self._debug(" Adding MOSEK conic constraint: {} in {}".format(
437 mosekVars, "Qr" if isRotated else "Q"))
439 # Add the conic constraint.
440 coneIndex = self.int.getnumcone()
441 mosekCone = mosek.conetype.rquad if isRotated else mosek.conetype.quad
442 self.int.appendcone(mosekCone, 0.0, mosekVars)
444 self._mosekCone[picosConstraint] = (coneIndex, mosekVars, mosekCons)
446 def _import_sdp_constraint(self, picosConstraint):
447 import mosek
449 n = picosConstraint.size[0]
450 dim = (n*(n + 1)) // 2
452 # MOSEK does not support general LMIs but so called "bar vars" which
453 # are variables in the symmetric positive semidefinite cone. We use them
454 # in combination with linear equalities to represent the LMI.
455 barVar = self.int.getnumbarvar()
456 mosekConOffset = self.int.getnumcon()
457 self.int.appendbarvars([n])
458 self.int.appendcons(dim)
460 # MOSEK uses a storage of symmetric coefficient matrices that are used
461 # as dot product coefficients to build scalar constraints involving both
462 # "bar vars" and normal scalar variables. We build a couple of these
463 # matrices to be able to select individual entries of our "bar vars".
464 # More precisely, if X is a n×n symmetric matrix and 0 ≤ k ≤ n*(n+1)/2,
465 # then <Units[n][k],X> = h(X)[k] where h refers to the lower-triangular,
466 # column-major half-vectorization of X.
467 if n in self._mosekBarUnitCoefs:
468 Units = self._mosekBarUnitCoefs[n]
469 else:
470 Units = self._mosekBarUnitCoefs[n] = [
471 self.int.appendsparsesymmat(
472 n, [row], [col], [1.0 if row == col else 0.5])
473 for row, col in self._low_tri_indices(n)]
475 # We iterate over the lower triangular scalar sub-expressions of the
476 # expression that the PICOS constraint states to be PSD, and constrain
477 # them to be eqal to the MOSEK "bar var" at the same index.
478 psd = picosConstraint.psd
479 psdRows = psd.sparse_rows(self._mosekVarOffset, lowerTriangle=True)
480 for lowTriIndex, (row, col) in enumerate(self._low_tri_indices(n)):
481 mosekConIndex = mosekConOffset + lowTriIndex
482 indices, coefficients, offset = psdRows[lowTriIndex]
484 # The lower-triangular entries in the PSD-constrained matrix …
485 self.int.putarow(mosekConIndex, indices, coefficients)
487 # … minus the corresponding bar var entries …
488 self.int.putbaraij(
489 mosekConIndex, barVar, [Units[lowTriIndex]], [-1.0])
491 # … should equal zero.
492 self.int.putconbound(
493 mosekConIndex, mosek.boundkey.fx, -offset, -offset)
495 if self._debug():
496 self._debug(" Index {} ({}, {}): indices = {}, coefs = {}"
497 .format(lowTriIndex, row, col, indices, coefficients))
499 self._mosekLMI[picosConstraint] = (barVar, mosekConOffset)
501 def _import_constraint(self, picosConstraint):
502 if self._debug():
503 self._debug("Importing Constraint: {}".format(picosConstraint))
505 if isinstance(picosConstraint, AffineConstraint):
506 self._import_linear_constraint(picosConstraint)
507 elif isinstance(picosConstraint, ConvexQuadraticConstraint):
508 self._import_quad_constraint(picosConstraint)
509 elif isinstance(picosConstraint, SOCConstraint) \
510 or isinstance(picosConstraint, RSOCConstraint):
511 self._import_quad_conic_constraint(picosConstraint)
512 elif isinstance(picosConstraint, LMIConstraint):
513 self._import_sdp_constraint(picosConstraint)
514 else:
515 assert isinstance(picosConstraint, DummyConstraint), \
516 "Unexpected constraint type: {}".format(
517 picosConstraint.__class__.__name__)
519 def _reset_objective(self):
520 # Reset affine part.
521 numVars = self.int.getnumvar()
522 self.int.putclist(range(numVars), [0.0]*numVars)
523 self.int.putcfix(0.0)
525 # Reset quadratic part.
526 self.int.putqobj([], [], [])
528 def _import_affine_objective(self, picosObjective):
529 mosekIndices = []
530 mosekCoefficients = []
532 for picosVar, picosCoef in picosObjective._linear_coefs.items():
533 for localIndex in range(picosVar.dim):
534 if picosCoef[localIndex]:
535 mosekIndex = self._mosekVarOffset[picosVar] + localIndex
536 mosekIndices.append(mosekIndex)
537 mosekCoefficients.append(picosCoef[localIndex])
539 self.int.putclist(mosekIndices, mosekCoefficients)
541 if picosObjective._constant_coef:
542 self.int.putcfix(picosObjective._constant_coef[0])
544 def _import_quadratic_objective(self, picosObjective):
545 # Import the quadratic part.
546 self.int.putqobj(*self._quadexp_pic2msk(picosObjective))
548 # Import the affine part.
549 self._import_affine_objective(picosObjective.aff)
551 def _import_objective(self):
552 import mosek
554 picosSense, picosObjective = self.ext.no
556 # Import objective sense.
557 if picosSense == "min":
558 self.int.putobjsense(mosek.objsense.minimize)
559 else:
560 assert picosSense == "max"
561 self.int.putobjsense(mosek.objsense.maximize)
563 # Import objective function.
564 if isinstance(picosObjective, AffineExpression):
565 self._import_affine_objective(picosObjective)
566 else:
567 assert isinstance(picosObjective, QuadraticExpression)
568 self._import_quadratic_objective(picosObjective)
570 def _import_problem(self):
571 # Create a problem instance.
572 self.int = self.env.Task()
574 # Import variables.
575 for variable in self.ext.variables.values():
576 self._import_variable(variable)
578 # Import constraints.
579 for constraint in self.ext.constraints.values():
580 self._import_constraint(constraint)
582 # Set objective.
583 self._import_objective()
585 def _update_problem(self):
586 for oldConstraint in self._removed_constraints():
587 raise ProblemUpdateError("PICOS does not support removing "
588 "constraints from a MOSEK instance.")
590 for oldVariable in self._removed_variables():
591 raise ProblemUpdateError("PICOS does not support removing variables"
592 " from a MOSEK instance.")
594 for newVariable in self._new_variables():
595 self._import_variable(newVariable)
597 for newConstraint in self._new_constraints():
598 self._import_constraint(newConstraint)
600 if self._objective_has_changed():
601 self._reset_objective()
602 self._import_objective()
604 def _solve(self):
605 import mosek
607 # Determine whether a basic solution will be available.
608 # TODO: Give Problem an interface for checks like this.
609 isLP = isinstance(self.ext.no.function, AffineExpression) \
610 and all(isinstance(constraint, AffineConstraint)
611 for constraint in self.ext.constraints.values())
613 # Reset all solver options to default.
614 self.int.setdefaults()
616 # verbosity
617 if self._verbose():
618 self.int.set_Stream(mosek.streamtype.log, self._streamprinter)
619 self.int.putintparam(mosek.iparam.log, self.ext.verbosity())
621 # abs_prim_fsb_tol
622 if self.ext.options.abs_prim_fsb_tol is not None:
623 value = self.ext.options.abs_prim_fsb_tol
625 # Interior-point primal feasibility tolerances.
626 self.int.putdouparam(mosek.dparam.intpnt_tol_pfeas, value)
627 self.int.putdouparam(mosek.dparam.intpnt_co_tol_pfeas, value)
628 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_pfeas, value)
629 if self.ver <= 8:
630 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_pfeas, value)
632 # Simplex primal feasibility tolerance.
633 self.int.putdouparam(mosek.dparam.basis_tol_x, value)
635 # Mixed-integer (primal) feasibility tolerance.
636 self.int.putdouparam(mosek.dparam.mio_tol_feas, value)
638 # abs_dual_fsb_tol
639 if self.ext.options.abs_dual_fsb_tol is not None:
640 value = self.ext.options.abs_dual_fsb_tol
642 # Interior-point dual feasibility tolerances.
643 self.int.putdouparam(mosek.dparam.intpnt_tol_dfeas, value)
644 self.int.putdouparam(mosek.dparam.intpnt_co_tol_dfeas, value)
645 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_dfeas, value)
646 if self.ver <= 8:
647 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_dfeas, value)
649 # Simplex dual feasibility (optimality) tolerance.
650 self.int.putdouparam(mosek.dparam.basis_tol_s, value)
652 # rel_dual_fsb_tol
653 if self.ext.options.rel_dual_fsb_tol is not None:
654 # Simplex relative dual feasibility (optimality) tolerance.
655 self.int.putdouparam(mosek.dparam.basis_rel_tol_s,
656 self.ext.options.rel_dual_fsb_tol)
658 # rel_ipm_opt_tol
659 if self.ext.options.rel_ipm_opt_tol is not None:
660 value = self.ext.options.rel_ipm_opt_tol
662 # Interior-point primal feasibility tolerances.
663 self.int.putdouparam(mosek.dparam.intpnt_tol_rel_gap, value)
664 self.int.putdouparam(mosek.dparam.intpnt_co_tol_rel_gap, value)
665 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_rel_gap, value)
666 if self.ver <= 8:
667 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_rel_gap, value)
669 # abs_bnb_opt_tol
670 if self.ext.options.abs_bnb_opt_tol is not None:
671 self.int.putdouparam(mosek.dparam.mio_tol_abs_gap,
672 self.ext.options.abs_bnb_opt_tol)
674 # rel_bnb_opt_tol
675 if self.ext.options.rel_bnb_opt_tol is not None:
676 self.int.putdouparam(mosek.dparam.mio_tol_rel_gap,
677 self.ext.options.rel_bnb_opt_tol)
679 # integrality_tol
680 if self.ext.options.integrality_tol is not None:
681 self.int.putdouparam(mosek.dparam.mio_tol_abs_relax_int,
682 self.ext.options.integrality_tol)
684 # max_iterations
685 if self.ext.options.max_iterations is not None:
686 value = self.ext.options.max_iterations
687 self.int.putintparam(mosek.iparam.bi_max_iterations, value)
688 self.int.putintparam(mosek.iparam.intpnt_max_iterations, value)
689 self.int.putintparam(mosek.iparam.sim_max_iterations, value)
691 _lpm = {
692 "interior": (mosek.optimizertype.intpnt if isLP
693 else mosek.optimizertype.conic),
694 "psimplex": mosek.optimizertype.primal_simplex,
695 "dsimplex": mosek.optimizertype.dual_simplex}
697 # lp_node_method
698 if self.ext.options.lp_node_method is not None:
699 value = self.ext.options.lp_node_method
700 assert value in _lpm, "Unexpected lp_node_method value."
701 self.int.putintparam(mosek.iparam.mio_node_optimizer, _lpm[value])
703 # lp_root_method
704 if self.ext.options.lp_root_method is not None:
705 value = self.ext.options.lp_root_method
706 assert value in _lpm, "Unexpected lp_root_method value."
707 self.int.putintparam(mosek.iparam.mio_root_optimizer, _lpm[value])
709 # timelimit
710 if self.ext.options.timelimit is not None:
711 value = float(self.ext.options.timelimit)
712 self.int.putdouparam(mosek.dparam.optimizer_max_time, value)
713 self.int.putdouparam(mosek.dparam.mio_max_time, value)
715 # max_fsb_nodes
716 if self.ext.options.max_fsb_nodes is not None:
717 self.int.putintparam(mosek.iparam.mio_max_num_solutions,
718 self.ext.options.max_fsb_nodes)
720 # Handle MOSEK-specific options.
721 for key, value in self.ext.options.mosek_params.items():
722 try:
723 self.int.putparam(key.upper(), str(value))
724 except mosek.Error as error:
725 self._handle_bad_solver_specific_option(key, value, error)
727 # Handle unsupported options.
728 # TODO: Handle "hotstart" option (via mio_construct_sol).
729 self._handle_unsupported_option("hotstart", "treememory")
731 # Attempt to solve the problem.
732 with self._header(), self._stopwatch():
733 try:
734 self.int.optimize()
735 except mosek.Error as error:
736 if error.errno in (
737 mosek.rescode.err_con_q_not_psd,
738 mosek.rescode.err_con_q_not_nsd,
739 mosek.rescode.err_obj_q_not_psd,
740 mosek.rescode.err_obj_q_not_nsd,
741 mosek.rescode.err_toconic_constr_q_not_psd,
742 mosek.rescode.err_toconic_objective_not_psd):
743 self._handle_continuous_nonconvex_error(error)
744 else:
745 raise
747 # Set the solution to be retrieved.
748 if self.ext.is_continuous():
749 if isLP:
750 solType = mosek.soltype.bas
751 else:
752 solType = mosek.soltype.itr
753 else:
754 solType = mosek.soltype.itg
756 # Retrieve primals.
757 primals = {}
758 if self.ext.options.primals is not False:
759 values = [float("nan")]*self.int.getnumvar()
760 self.int.getxx(solType, values)
762 for picosVar in self.ext.variables.values():
763 mosekOffset = self._mosekVarOffset[picosVar]
764 primal = values[mosekOffset:mosekOffset + picosVar.dim]
766 if float("nan") in primal:
767 primals[picosVar] = None
768 else:
769 primals[picosVar] = primal
771 # Retrieve duals.
772 duals = {}
773 if self.ext.options.duals is not False and self.ext.is_continuous():
774 minimizing = self.int.getobjsense() == mosek.objsense.minimize
776 for constraint in self.ext.constraints.values():
777 if isinstance(constraint, DummyConstraint):
778 duals[constraint] = cvxopt.spmatrix(
779 [], [], [], constraint.size)
780 continue
782 length = len(constraint)
783 dual = [float("nan")]*length
785 if isinstance(constraint, AffineConstraint):
786 offset = self._mosekLinConOffset[constraint]
787 self.int.getyslice(solType, offset, offset + length, dual)
788 elif isinstance(constraint, ConvexQuadraticConstraint):
789 # TODO: Implement consistent QCQP dual retrieval for all
790 # solvers that return such duals.
791 dual = None
792 elif isinstance(constraint, SOCConstraint) \
793 or isinstance(constraint, RSOCConstraint):
794 mosekVars = self._mosekCone[constraint][1]
795 for localConeIndex in range(length):
796 x = [float("nan")]
797 offset = mosekVars[localConeIndex]
798 self.int.getsnxslice(solType, offset, offset + 1, x)
799 dual[localConeIndex] = x[0]
801 if isinstance(constraint, SOCConstraint):
802 dual = [-du for du in dual]
803 elif isinstance(constraint, RSOCConstraint):
804 dual[0] = -dual[0] / math.sqrt(2.0)
805 dual[1] = -dual[1] / math.sqrt(2.0)
806 dual[2:] = [-du for du in dual[2:]]
807 elif isinstance(constraint, LMIConstraint):
808 n = constraint.size[0]
809 barVar, _ = self._mosekLMI[constraint]
810 lowerTriangularDual = [float("nan")]*(n*(n + 1) // 2)
811 self.int.getbarsj(solType, barVar, lowerTriangularDual)
812 for lti, (row, col) in enumerate(self._low_tri_indices(n)):
813 value = lowerTriangularDual[lti]
814 dual[n*row + col] = value
815 if row != col:
816 dual[n*col + row] = value
817 else:
818 assert False, "Constraint type not supported."
820 if dual is None:
821 pass
822 elif float("nan") in dual:
823 dual = None
824 else:
825 dual = cvxopt.matrix(dual, constraint.size)
827 if isinstance(constraint, AffineConstraint) \
828 or isinstance(constraint, LMIConstraint):
829 if not constraint.is_increasing():
830 dual = -dual
832 if minimizing:
833 dual = -dual
835 duals[constraint] = dual
837 # Retrieve objective value.
838 value = self.int.getprimalobj(solType)
840 # Retrieve solution status.
841 primalStatus = self._solution_status_pic2msk(
842 self.int.getsolsta(solType), primalSolution=True)
843 dualStatus = self._solution_status_pic2msk(
844 self.int.getsolsta(solType), primalSolution=False)
845 problemStatus = self._problem_status_pic2msk(
846 self.int.getprosta(solType), not self.ext.is_continuous())
848 return self._make_solution(
849 value, primals, duals, primalStatus, dualStatus, problemStatus)
851 def _solution_status_pic2msk(self, statusCode, primalSolution):
852 from mosek import solsta as ss
853 dualSolution = not primalSolution
855 map = {
856 ss.unknown: SS_UNKNOWN,
857 ss.optimal: SS_OPTIMAL,
858 ss.prim_feas:
859 SS_FEASIBLE if primalSolution else SS_UNKNOWN,
860 ss.dual_feas:
861 SS_FEASIBLE if dualSolution else SS_UNKNOWN,
862 ss.prim_and_dual_feas: SS_FEASIBLE,
863 ss.prim_infeas_cer:
864 SS_INFEASIBLE if primalSolution else SS_UNKNOWN,
865 ss.dual_infeas_cer:
866 SS_INFEASIBLE if dualSolution else SS_UNKNOWN,
867 ss.prim_illposed_cer: SS_UNKNOWN,
868 ss.dual_illposed_cer: SS_UNKNOWN,
869 ss.integer_optimal: SS_OPTIMAL,
870 }
872 if self.ver < 9:
873 map.update({
874 ss.near_optimal: SS_FEASIBLE,
875 ss.near_prim_feas: SS_UNKNOWN,
876 ss.near_dual_feas: SS_UNKNOWN,
877 ss.near_prim_and_dual_feas: SS_UNKNOWN,
878 ss.near_prim_infeas_cer: SS_UNKNOWN,
879 ss.near_dual_infeas_cer: SS_UNKNOWN,
880 ss.near_integer_optimal: SS_FEASIBLE
881 })
883 try:
884 return map[statusCode]
885 except KeyError:
886 self._warn(
887 "The MOSEK solution status code {} is not known to PICOS."
888 .format(statusCode))
889 return SS_UNKNOWN
891 def _problem_status_pic2msk(self, statusCode, integerProblem):
892 from mosek import prosta as ps
894 map = {
895 ps.unknown: PS_UNKNOWN,
896 ps.prim_and_dual_feas: PS_FEASIBLE,
897 ps.prim_feas:
898 PS_FEASIBLE if integerProblem else PS_UNKNOWN,
899 ps.dual_feas: PS_UNKNOWN,
900 ps.prim_infeas: PS_INFEASIBLE,
901 ps.dual_infeas: PS_UNBOUNDED, # TODO: UNB_OR_INF
902 ps.prim_and_dual_infeas: PS_INFEASIBLE,
903 ps.ill_posed: PS_ILLPOSED,
904 ps.prim_infeas_or_unbounded: PS_UNKNOWN # TODO: UNB_OR_INF
905 }
907 if self.ver < 9:
908 map.update({
909 ps.near_prim_and_dual_feas: PS_UNKNOWN,
910 ps.near_prim_feas: PS_UNKNOWN,
911 ps.near_dual_feas: PS_UNKNOWN
912 })
914 try:
915 return map[statusCode]
916 except KeyError:
917 self._warn("The MOSEK problem status code {} is not known to PICOS."
918 .format(statusCode))
919 return PS_UNKNOWN
922# --------------------------------------
923__all__ = api_end(_API_START, globals())