Hide keyboard shortcuts

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# ------------------------------------------------------------------------------ 

18 

19"""Implementation of :class:`MOSEKSolver`.""" 

20 

21import math 

22import sys 

23 

24import cvxopt 

25 

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 

37 

38_API_START = api_start(globals()) 

39# ------------------------------- 

40 

41 

42class MOSEKSolver(Solver): 

43 """Interface to the MOSEK solver via its low level Optimizer API. 

44 

45 Supports both MOSEK 8 and 9. 

46 

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 """ 

50 

51 SUPPORTED = Specification( 

52 objectives=[ 

53 AffineExpression, 

54 QuadraticExpression], 

55 constraints=[ 

56 DummyConstraint, 

57 AffineConstraint, 

58 SOCConstraint, 

59 RSOCConstraint, 

60 ConvexQuadraticConstraint, 

61 LMIConstraint]) 

62 

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 

69 

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 

76 

77 conic = any(("con", constraint) in footprint 

78 for constraint in (SOCConstraint, RSOCConstraint, LMIConstraint)) 

79 

80 quadratic = footprint.objective.clstype is QuadraticExpression \ 

81 or ("con", ConvexQuadraticConstraint) in footprint 

82 

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 

90 

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 

97 

98 if footprint not in cls.SUPPORTED: 

99 if explain: 

100 return False, cls.SUPPORTED.mismatch_reason(footprint) 

101 else: 

102 return False 

103 

104 return (True, None) if explain else True 

105 

106 @classmethod 

107 def default_penalty(cls): 

108 """Implement :meth:`~.solver.Solver.default_penalty`.""" 

109 return 0.0 # Commercial solver. 

110 

111 @classmethod 

112 def test_availability(cls): 

113 """Implement :meth:`~.solver.Solver.test_availability`.""" 

114 cls.check_import("mosek") 

115 

116 @classmethod 

117 def names(cls): 

118 """Implement :meth:`~.solver.Solver.names`.""" 

119 return "mosek", "MOSEK (Optimizer)", "MOSEK via Optimizer API" 

120 

121 @classmethod 

122 def is_free(cls): 

123 """Implement :meth:`~.solver.Solver.is_free`.""" 

124 return False 

125 

126 def __init__(self, problem): 

127 """Initialize a MOSEK (Optimizer) solver interface. 

128 

129 :param ~picos.Problem problem: The problem to be solved. 

130 """ 

131 super(MOSEKSolver, self).__init__(problem) 

132 

133 self._mosekVarOffset = dict() 

134 """Maps a PICOS (multidimensional) variable to a MOSEK scalar variable 

135 (start) index.""" 

136 

137 self._mosekBarVarIndex = dict() 

138 """Maps suited PICOS symmetric variables to a MOSEK bar var index.""" 

139 

140 self._mosekLinConOffset = dict() 

141 """Maps a PICOS (multidimensional) linear constraint to a MOSEK scalar 

142 constraint (start) index.""" 

143 

144 self._mosekQuadConIndex = dict() 

145 """Maps a PICOS quadratic constraint to a MOSEK constraint index.""" 

146 

147 self._mosekCone = dict() 

148 """ 

149 Maps a PICOS SOCC or RSOCC to a triple: 

150 

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 """ 

162 

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.""" 

167 

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.""" 

175 

176 def reset_problem(self): 

177 """Implement :meth:`~.solver.Solver.reset_problem`.""" 

178 self.int = None 

179 

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() 

187 

188 @classmethod 

189 def _get_environment(cls): 

190 if not hasattr(cls, "mosekEnvironment"): 

191 import mosek 

192 cls.mosekEnvironment = mosek.Env() 

193 

194 return cls.mosekEnvironment 

195 

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.")""" 

200 

201 @classmethod 

202 def _get_major_version(cls): 

203 if not hasattr(cls, "mosekVersion"): 

204 import mosek 

205 cls.mosekVersion = mosek.Env.getversion() 

206 

207 return cls.mosekVersion[0] 

208 

209 ver = property(lambda self: self.__class__._get_major_version()) 

210 """The major version of the available MOSEK library.""" 

211 

212 @staticmethod 

213 def _streamprinter(text): 

214 sys.stdout.write(text) 

215 sys.stdout.flush() 

216 

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 

225 

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 

233 

234 def _affinexp_pic2msk(self, picosExpression): 

235 """Convert a PICOS affine expression to a MOSEK one. 

236 

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 

249 

250 return 

251 

252 barVarIndices = tuple([] for _ in range(len(picosExpression))) 

253 barVarCoefs = tuple([] for _ in range(len(picosExpression))) 

254 

255 for var, coef in picosExpression._linear_coefs.items(): 

256 if var in self._mosekBarVarIndex: 

257 mosekBarVarIndex = self._mosekBarVarIndex[var] 

258 

259 assert isinstance(var, SymmetricVariable) 

260 

261 n = var.shape[0] 

262 

263 for localIndex in range(coef.size[0]): 

264 localCoef = coef[localIndex, :] 

265 

266 if not localCoef: 

267 continue 

268 

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) 

273 

274 assert localCoefMatrix.size == (n, n) 

275 

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) 

286 

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) 

292 

293 # Append bar var index and coefficient matrix index. 

294 barVarIndices[localIndex].append(mosekBarVarIndex) 

295 barVarCoefs[localIndex].append(mosekCoefMatrix) 

296 

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))) 

301 

302 J, V = JV if JV else ((),)*2 

303 K, W = barVarIndices[localIndex], barVarCoefs[localIndex] 

304 

305 yield J, V, K, W, c 

306 

307 def _scalar_affinexp_pic2msk(self, picosExpression): 

308 assert len(picosExpression) == 1 

309 return next(self._affinexp_pic2msk(picosExpression)) 

310 

311 def _quadexp_pic2msk(self, picosExpression): 

312 """Transform a quadratic expression from PICOS to MOSEK. 

313 

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() 

322 

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 

326 

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 

340 

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 

347 

348 I.append(mskVar1Index) 

349 J.append(mskVar2Index) 

350 V.append(localCoefficient) 

351 

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 

357 

358 # Return it as a sparse triple for MOSEK to consume. 

359 return list(Q.I), list(Q.J), list(Q.V) 

360 

361 def _import_variable(self, picosVar): 

362 import mosek 

363 

364 numVars = self._mosekVarOffset[picosVar] = self.int.getnumvar() 

365 dim = picosVar.dim 

366 indices = range(numVars, numVars + dim) 

367 self.int.appendvars(dim) 

368 

369 # Set the variable type. 

370 if isinstance(picosVar, (IntegerVariable, BinaryVariable)): 

371 self.int.putvartypelist( 

372 indices, [mosek.variabletype.type_int]*dim) 

373 

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 

382 

383 lower, upper = picosVar.bound_dicts 

384 

385 for i, b in lower.items(): 

386 boundKeys[i] = mosek.boundkey.lo 

387 lowerBounds[i] = b 

388 

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 

397 

398 upperBounds[i] = b 

399 

400 self.int.putvarboundlist( 

401 indices, boundKeys, lowerBounds, upperBounds) 

402 

403 def _import_psd_variable(self, picosVar): 

404 assert isinstance(picosVar, SymmetricVariable) 

405 

406 self._mosekBarVarIndex[picosVar] = self.int.getnumbarvar() 

407 self.int.appendbarvars([picosVar.shape[0]]) 

408 

409 def _import_linear_constraint(self, picosConstraint): 

410 import mosek 

411 

412 numCons = self.int.getnumcon() 

413 conLen = len(picosConstraint) 

414 self.int.appendcons(conLen) 

415 

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 

422 

423 rows = self._affinexp_pic2msk(picosConstraint.lhs - picosConstraint.rhs) 

424 

425 for localConIndex, (J, V, K, W, c) in enumerate(rows): 

426 mosekConIndex = numCons + localConIndex 

427 rhs = -c 

428 

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) 

433 

434 self._mosekLinConOffset[picosConstraint] = numCons 

435 

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 

440 

441 # Import the linear part first. 

442 picosLinConPart = picosConstraint.le0.aff < 0 

443 self._import_linear_constraint(picosLinConPart) 

444 mosekConIndex = self._mosekLinConOffset.pop(picosLinConPart) 

445 

446 # Add the quadratic part. 

447 self.int.putqconk( 

448 mosekConIndex, *self._quadexp_pic2msk(picosConstraint.le0)) 

449 

450 self._mosekQuadConIndex[picosConstraint] = mosekConIndex 

451 

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 

459 

460 def _import_quad_conic_constraint(self, picosConstraint): 

461 import mosek 

462 

463 isRotated = isinstance(picosConstraint, RSOCConstraint) 

464 mosekVars, mosekCons = [], [None]*len(picosConstraint) 

465 

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)) 

478 

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) 

488 

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) 

495 

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) 

499 

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 

507 

508 # Prepare the auxiliary constraint. 

509 J, V, K, W, c = entries[missingVarIndex] 

510 

511 if self._debug(): 

512 self._debug(" Adding MOSEK auxiliary constraint: " 

513 "{}.T * x{}{} = {}" 

514 .format(V, J, " + f(barvars)" if K else "", -c)) 

515 

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) 

522 

523 # Complete the mapping of cone member entries to MOSEK (auxiliary) 

524 # variables (and auxiliary constraints). 

525 mosekVars[missingVarIndex] = auxVarIndex 

526 mosekCons[missingVarIndex] = auxConIndex 

527 

528 if self._debug(): 

529 self._debug(" Adding MOSEK conic constraint: {} in {}".format( 

530 mosekVars, "Qr" if isRotated else "Q")) 

531 

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) 

536 

537 self._mosekCone[picosConstraint] = (coneIndex, mosekVars, mosekCons) 

538 

539 def _import_sdp_constraint(self, picosConstraint): 

540 # NOTE: Trivial LMIs of the form X ≽ 0 are loaded as bar vars instead. 

541 

542 import mosek 

543 

544 n = picosConstraint.size[0] 

545 dim = (n*(n + 1)) // 2 

546 

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) 

554 

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)] 

569 

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 

578 

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]) 

583 

584 # … minus the corresponding bar var entries … 

585 self.int.putbaraij( 

586 mosekConIndex, barVar, [Units[lowTriIndex]], [-1.0]) 

587 

588 # … should equal zero. 

589 self.int.putconbound(mosekConIndex, mosek.boundkey.fx, rhs, rhs) 

590 

591 if self._debug(): 

592 self._debug(" Index {} ({}, {}): J = {}, V = {}, ..." 

593 .format(lowTriIndex, row, col, J, V)) 

594 

595 self._mosekLMI[picosConstraint] = (barVar, mosekConOffset) 

596 

597 def _import_constraint(self, picosConstraint): 

598 if self._debug(): 

599 self._debug("Importing Constraint: {}".format(picosConstraint)) 

600 

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__) 

614 

615 def _reset_objective(self): 

616 numVars = self.int.getnumvar() 

617 numBarVars = self.int.getnumbarvar() 

618 

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) 

624 

625 # Reset quadratic part. 

626 self.int.putqobj([], [], []) 

627 

628 def _import_affine_objective(self, picosObjective): 

629 J, V, K, W, c = self._scalar_affinexp_pic2msk(picosObjective) 

630 

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) 

635 

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 

640 

641 # Import the quadratic part. 

642 self.int.putqobj(*self._quadexp_pic2msk(picosObjective)) 

643 

644 # Import the affine part. 

645 self._import_affine_objective(picosObjective.aff) 

646 

647 def _import_objective(self): 

648 import mosek 

649 

650 picosSense, picosObjective = self.ext.no 

651 

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) 

658 

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) 

665 

666 def _import_problem(self): 

667 # Create a problem instance. 

668 self.int = self.env.Task() 

669 

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) 

678 

679 # Import normal variables. 

680 for variable in self.ext.variables.values(): 

681 if variable not in psdVars: 

682 self._import_variable(variable) 

683 

684 # Import virtual PSD variables as MOSEK "bar variables". 

685 for variable in psdVars: 

686 self._import_psd_variable(variable) 

687 

688 # Import remaining constraints. 

689 for constraint in constraints: 

690 self._import_constraint(constraint) 

691 

692 # Set objective. 

693 self._import_objective() 

694 

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.") 

699 

700 for oldVariable in self._removed_variables(): 

701 raise ProblemUpdateError("PICOS does not support removing variables" 

702 " from a MOSEK instance.") 

703 

704 for newVariable in self._new_variables(): 

705 self._import_variable(newVariable) 

706 

707 for newConstraint in self._new_constraints(): 

708 self._import_constraint(newConstraint) 

709 

710 if self._objective_has_changed(): 

711 self._reset_objective() 

712 self._import_objective() 

713 

714 def _solve(self): 

715 import mosek 

716 

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()) 

722 

723 # Reset all solver options to default. 

724 self.int.setdefaults() 

725 self.int.putoptserverhost("") 

726 

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()) 

731 

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 

735 

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) 

742 

743 # Simplex primal feasibility tolerance. 

744 self.int.putdouparam(mosek.dparam.basis_tol_x, value) 

745 

746 # Mixed-integer (primal) feasibility tolerance. 

747 self.int.putdouparam(mosek.dparam.mio_tol_feas, value) 

748 

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 

752 

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) 

759 

760 # Simplex dual feasibility (optimality) tolerance. 

761 self.int.putdouparam(mosek.dparam.basis_tol_s, value) 

762 

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) 

768 

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 

772 

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) 

779 

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) 

784 

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) 

789 

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) 

794 

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) 

801 

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} 

807 

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]) 

813 

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]) 

819 

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) 

825 

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) 

830 

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) 

837 

838 # Handle 'mosek_server' option. 

839 if self.ext.options.mosek_server: 

840 self.int.putoptserverhost(self.ext.options.mosek_server) 

841 

842 # Handle unsupported options. 

843 # TODO: Handle "hotstart" option (via mio_construct_sol). 

844 self._handle_unsupported_option("hotstart", "treememory") 

845 

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 

861 

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 

870 

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) 

876 

877 for picosVar in self.ext.variables.values(): 

878 if picosVar in self._mosekBarVarIndex: 

879 barVarIndex = self._mosekBarVarIndex[picosVar] 

880 

881 assert isinstance(picosVar, SymmetricVariable) 

882 

883 n, d = picosVar.shape[0], picosVar.dim 

884 lowTriPrimal = [float("nan")]*d 

885 matrixPrimal = [float("nan")]*n**2 

886 

887 # Obtain a lower triangular vectorization from MOSEK. 

888 self.int.getbarxj(solType, barVarIndex, lowTriPrimal) 

889 

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 

897 

898 # Load the full vectorization as a CVXOPT matrix. 

899 matrixPrimal = cvxopt.matrix(matrixPrimal, (n, n)) 

900 

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] 

906 

907 if float("nan") in primal: 

908 primals[picosVar] = None 

909 else: 

910 primals[picosVar] = primal 

911 

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 

916 

917 for constraint in self.ext.constraints.values(): 

918 if isinstance(constraint, DummyConstraint): 

919 duals[constraint] = cvxopt.spmatrix( 

920 [], [], [], constraint.size) 

921 continue 

922 

923 length = len(constraint) 

924 dual = [float("nan")]*length 

925 

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] 

941 

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] 

958 

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." 

970 

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) 

977 

978 if isinstance(constraint, AffineConstraint) \ 

979 or isinstance(constraint, LMIConstraint): 

980 if not constraint.is_increasing(): 

981 dual = -dual 

982 

983 if minimizing: 

984 dual = -dual 

985 

986 duals[constraint] = dual 

987 

988 # Retrieve objective value. 

989 value = self.int.getprimalobj(solType) 

990 

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()) 

998 

999 return self._make_solution( 

1000 value, primals, duals, primalStatus, dualStatus, problemStatus) 

1001 

1002 def _solution_status_pic2msk(self, statusCode, primalSolution): 

1003 from mosek import solsta as ss 

1004 dualSolution = not primalSolution 

1005 

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 } 

1022 

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 }) 

1033 

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 

1041 

1042 def _problem_status_pic2msk(self, statusCode, integerProblem): 

1043 from mosek import prosta as ps 

1044 

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 } 

1057 

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 }) 

1064 

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 

1071 

1072 

1073# -------------------------------------- 

1074__all__ = api_end(_API_START, globals())