Coverage for picos/solvers/solver_mosek.py: 88.61%

553 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-03-26 07:46 +0000

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", "MOSEK", "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 variable to a MOSEK scalar variable offset.""" 

135 

136 self._mosekBarVarIndex = dict() 

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

138 

139 self._mosekLinConOffset = dict() 

140 """Maps a PICOS linear constraint to a MOSEK constraint offset.""" 

141 

142 self._mosekQuadConIndex = dict() 

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

144 

145 self._mosekCone = dict() 

146 """ 

147 Maps a PICOS SOCC or RSOCC to a triple: 

148 

149 - The first entry is the index of a MOSEK conic constraint. 

150 - The second entry is a list of MOSEK scalar (auxiliary) variable 

151 indices that appear in the MOSEK conic constraint. 

152 - The third entry is another list of same size containing auxiliary 

153 constraints (or None). If an entry in the second list is None, then 

154 the corresponding index in the first list is of a proper MOSEK scalar 

155 variable instead of an auxiliary one, which can happen at most once 

156 per variable and only if the structure of the PICOS constraint allows 

157 it (that is, if the repspective entry of the cone member happens to be 

158 a PICOS scalar variable as opposed to a composed affine expression). 

159 """ 

160 

161 self._mosekLMI = dict() 

162 """Maps a PICOS LMI to a pair representing its MOSEK representation. The 

163 first entry is the index of a MOSEK symmetric PSD "bar variable", the 

164 second entry is the offset for a number of scalar linear equalities.""" 

165 

166 self._mosekBarUnitCoefs = dict() 

167 """Maps the row (column) count of a symmetric matrix to a list of MOSEK 

168 symmetric coefficient matrices. More precisely, if X is a n×n symmetric 

169 matrix, 0 ≤ k ≤ n*(n+1)/2, then <_mosekBarUnitCoefs[n][k], X> = h(X)[k] 

170 where h refers to the lower-triangular, column-major half-vectorization 

171 of X. These matrices are used to represent LMIs. Unlike values of 

172 _mosekLMI, they are shared between LMIs of same size.""" 

173 

174 def reset_problem(self): 

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

176 self.int = None 

177 

178 self._mosekVarOffset.clear() 

179 self._mosekBarVarIndex.clear() 

180 self._mosekLinConOffset.clear() 

181 self._mosekQuadConIndex.clear() 

182 self._mosekCone.clear() 

183 self._mosekLMI.clear() 

184 self._mosekBarUnitCoefs.clear() 

185 

186 @classmethod 

187 def _get_environment(cls): 

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

189 import mosek 

190 cls.mosekEnvironment = mosek.Env() 

191 

192 return cls.mosekEnvironment 

193 

194 env = property(lambda self: self.__class__._get_environment()) 

195 """This references a MOSEK environment, which is shared among all 

196 MOSEKSolver instances. (The MOSEK documentation states that "[a]ll tasks in 

197 the program should share the same environment.")""" 

198 

199 @classmethod 

200 def _get_major_version(cls): 

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

202 import mosek 

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

204 

205 return cls.mosekVersion[0] 

206 

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

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

209 

210 @staticmethod 

211 def _streamprinter(text): 

212 sys.stdout.write(text) 

213 sys.stdout.flush() 

214 

215 @staticmethod 

216 def _low_tri_indices(rowCount): 

217 """Yield enumerated lower triangular indices in column-major order.""" 

218 lti = 0 

219 for col in range(rowCount): 

220 for row in range(col, rowCount): 

221 yield lti, row, col 

222 lti += 1 

223 

224 def _affinexp_pic2msk(self, picosExpression): 

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

226 

227 :yield: 

228 A sequence of quintuples ``(J, V, K, W, c)``, each corresponding to 

229 one scalar entry of a multidimensional PICOS affine expression, 

230 where ``J`` are MOSEK scalar variable indices, ``V`` are associated 

231 coefficients, ``K`` are MOSEK bar var (symmetric positive 

232 semidefinite variable) indices, ``W`` are indices of associated 

233 coefficient matrices and where ``c`` is a scalar constant. 

234 """ 

235 # Handle the base case (no bar vars in the problem) quickly. 

236 if not self._mosekBarVarIndex: 

237 for J, V, c in picosExpression.sparse_rows(self._mosekVarOffset): 

238 yield J, V, (), (), c 

239 

240 return 

241 

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

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

244 

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

246 if var in self._mosekBarVarIndex: 

247 mosekBarVarIndex = self._mosekBarVarIndex[var] 

248 

249 assert isinstance(var, SymmetricVariable) 

250 

251 n = var.shape[0] 

252 

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

254 localCoef = coef[localIndex, :] 

255 

256 if not localCoef: 

257 continue 

258 

259 # Devectorize local coefficient. 

260 # This works as c.T * svec(X) = <desvec(c), X> due to svec 

261 # and its inverse desvec being isometric isomorphisms. 

262 localCoefMatrix = var._vec.devectorize(localCoef.T) 

263 

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

265 

266 # Get lower triangular sparse rows for the coefficient. 

267 # NOTE: This is much faster than using _low_low_tri_indices. 

268 I, J, V = [], [], [] 

269 for j in range(n): 

270 for i in range(j, n): 

271 v = localCoefMatrix[i, j] 

272 if v: 

273 I.append(i) 

274 J.append(j) 

275 V.append(v) 

276 

277 # Create a new MOSEK symmetric coefficient matrix. 

278 # NOTE: It does not seem possible to remove these matrices 

279 # in MOSEK 9.2 so we might leak memory whenever they 

280 # fall out of use (e.g. objective function updates). 

281 mosekCoefMatrix = self.int.appendsparsesymmat(n, I, J, V) 

282 

283 # Append bar var index and coefficient matrix index. 

284 barVarIndices[localIndex].append(mosekBarVarIndex) 

285 barVarCoefs[localIndex].append(mosekCoefMatrix) 

286 

287 rows = enumerate(picosExpression.sparse_rows(self._mosekVarOffset)) 

288 for localIndex, (J, V, c) in rows: 

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

290 

291 yield J, V, K, W, c 

292 

293 def _scalar_affinexp_pic2msk(self, picosExpression): 

294 assert len(picosExpression) == 1 

295 return next(self._affinexp_pic2msk(picosExpression)) 

296 

297 def _quadexp_pic2msk(self, picosExpression): 

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

299 

300 Tranforms the quadratic part of a PICOS quadratic expression to a 

301 symmetric, sparse biliniar form of which only the lower triangular 

302 entries are given, and that can be used with MOSEK's variable vector. 

303 Note that MOSEK applies a built-in factor of 0.5 to all biliniar forms 

304 while PICOS doesn't, so a factor of 2 is applied here to cancel it out. 

305 """ 

306 assert isinstance(picosExpression, QuadraticExpression) 

307 numVars = self.int.getnumvar() 

308 

309 # Mixing of conic and quadratic constraints is not supported by MOSEK; 

310 # therefor we do not handle bar vars here. 

311 assert not self._mosekBarVarIndex 

312 

313 # Make a sparse representation of the strict lower, diagonal and strict 

314 # upper parts of the matrix. 

315 IL, JL, VL = [], [], [] 

316 ID, JD, VD = [], [], [] 

317 IU, JU, VU = [], [], [] 

318 for (picosVar1, picosVar2), picosCoefs \ 

319 in picosExpression._sparse_quads.items(): 

320 for sparseIndex in range(len(picosCoefs)): 

321 localVar1Index = picosCoefs.I[sparseIndex] 

322 localVar2Index = picosCoefs.J[sparseIndex] 

323 localCoefficient = picosCoefs.V[sparseIndex] 

324 mskVar1Index = self._mosekVarOffset[picosVar1] + localVar1Index 

325 mskVar2Index = self._mosekVarOffset[picosVar2] + localVar2Index 

326 

327 if mskVar2Index < mskVar1Index: 

328 I, J, V = IL, JL, VL 

329 elif mskVar1Index < mskVar2Index: 

330 I, J, V = IU, JU, VU 

331 else: 

332 I, J, V = ID, JD, VD 

333 

334 I.append(mskVar1Index) 

335 J.append(mskVar2Index) 

336 V.append(localCoefficient) 

337 

338 # Compute the lower triangular part of the biliniar form. 

339 L = cvxopt.spmatrix(VL, IL, JL, (numVars, numVars)) 

340 D = cvxopt.spmatrix(VD, ID, JD, (numVars, numVars)) 

341 U = cvxopt.spmatrix(VU, IU, JU, (numVars, numVars)) 

342 Q = 2*D + L + U.T 

343 

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

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

346 

347 def _import_variable(self, picosVar): 

348 import mosek 

349 

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

351 dim = picosVar.dim 

352 indices = range(numVars, numVars + dim) 

353 self.int.appendvars(dim) 

354 

355 # Set the variable type. 

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

357 self.int.putvartypelist( 

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

359 

360 # Import bounds, including the implied bounds of a BinaryVariable. 

361 if isinstance(picosVar, BinaryVariable): 

362 self.int.putvarboundlist( 

363 indices, [mosek.boundkey.ra]*dim, [0]*dim, [1]*dim) 

364 else: 

365 boundKeys = [mosek.boundkey.fr]*dim 

366 lowerBounds = [0.0]*dim 

367 upperBounds = [0.0]*dim 

368 

369 lower, upper = picosVar.bound_dicts 

370 

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

372 boundKeys[i] = mosek.boundkey.lo 

373 lowerBounds[i] = b 

374 

375 for i, b in upper.items(): 

376 if boundKeys[i] == mosek.boundkey.fr: 

377 boundKeys[i] = mosek.boundkey.up 

378 else: # Also has a lower bound. 

379 if lowerBounds[i] == b: 

380 boundKeys[i] = mosek.boundkey.fx 

381 else: 

382 boundKeys[i] = mosek.boundkey.ra 

383 

384 upperBounds[i] = b 

385 

386 self.int.putvarboundlist( 

387 indices, boundKeys, lowerBounds, upperBounds) 

388 

389 def _import_psd_variable(self, picosVar): 

390 assert isinstance(picosVar, SymmetricVariable) 

391 

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

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

394 

395 def _import_linear_constraint(self, picosConstraint): 

396 import mosek 

397 

398 numCons = self.int.getnumcon() 

399 conLen = len(picosConstraint) 

400 self.int.appendcons(conLen) 

401 

402 if picosConstraint.is_equality(): 

403 boundKey = mosek.boundkey.fx 

404 elif picosConstraint.is_increasing(): 

405 boundKey = mosek.boundkey.up 

406 else: 

407 boundKey = mosek.boundkey.lo 

408 

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

410 

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

412 mosekConIndex = numCons + localConIndex 

413 rhs = -c 

414 

415 self.int.putarow(mosekConIndex, J, V) 

416 for k, w in zip(K, W): 

417 self.int.putbaraij(mosekConIndex, k, [w], [1.0]) 

418 self.int.putconbound(mosekConIndex, boundKey, rhs, rhs) 

419 

420 self._mosekLinConOffset[picosConstraint] = numCons 

421 

422 def _import_quad_constraint(self, picosConstraint): 

423 # Mixing of conic and quadratic constraints is not supported by MOSEK; 

424 # therefor we do not handle bar vars here. 

425 assert not self._mosekBarVarIndex 

426 

427 # Import the linear part first. 

428 picosLinConPart = picosConstraint.le0.aff < 0 

429 self._import_linear_constraint(picosLinConPart) 

430 mosekConIndex = self._mosekLinConOffset.pop(picosLinConPart) 

431 

432 # Add the quadratic part. 

433 self.int.putqconk( 

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

435 

436 self._mosekQuadConIndex[picosConstraint] = mosekConIndex 

437 

438 def _var_was_used_in_cone(self, mosekVariableIndex, usedJustNow=[]): 

439 if mosekVariableIndex in usedJustNow: 

440 return True 

441 for _, mosekVarIndices, _ in self._mosekCone.values(): 

442 if mosekVariableIndex in mosekVarIndices: 

443 return True 

444 return False 

445 

446 def _import_quad_conic_constraint(self, picosConstraint): 

447 import mosek 

448 

449 isRotated = isinstance(picosConstraint, RSOCConstraint) 

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

451 

452 # Get an initial MOSEK representation of the cone member. 

453 entries = [] 

454 if isRotated: 

455 # MOSEK internally adds a factor of 2 to the upper bound while PICOS 

456 # doesn't, so cancel it out by adding a factor of 1/sqrt(2) to both 

457 # factors of the upper bound. 

458 f = 1.0 / math.sqrt(2.0) 

459 entries.append(self._scalar_affinexp_pic2msk(f*picosConstraint.ub1)) 

460 entries.append(self._scalar_affinexp_pic2msk(f*picosConstraint.ub2)) 

461 else: 

462 entries.append(self._scalar_affinexp_pic2msk(picosConstraint.ub)) 

463 entries.extend(self._affinexp_pic2msk(picosConstraint.ne)) 

464 

465 # Map cone member entries to existing MOSEK variables, if possible. 

466 mosekVarsMissing = [] 

467 for scalarVarNum, (J, V, K, W, c) in enumerate(entries): 

468 if len(J) == 1 and V[0] == 1.0 and not K and not W and not c \ 

469 and not self._var_was_used_in_cone(J[0], mosekVars): 

470 mosekVars.append(J[0]) 

471 else: 

472 mosekVars.append(None) 

473 mosekVarsMissing.append(scalarVarNum) 

474 

475 # Create auxiliary variables and constraints. 

476 numAux = len(mosekVarsMissing) 

477 auxVarOffset = self.int.getnumvar() 

478 auxConOffset = self.int.getnumcon() 

479 self.int.appendvars(numAux) 

480 self.int.appendcons(numAux) 

481 

482 # Mosek fixes (!) new variables at zero, so set them free. 

483 self.int.putvarboundlist(range(auxVarOffset, auxVarOffset + numAux), 

484 [mosek.boundkey.fr]*numAux, [0.0]*numAux, [0.0]*numAux) 

485 

486 # Constrain the auxiliary variables to be equal to the cone member 

487 # entries for which no existing MOSEK variable could be used. 

488 # TODO: Instead of always creating a constraint, fix variables via their 

489 # bound whenever possible. 

490 for auxNum, missingVarIndex in enumerate(mosekVarsMissing): 

491 auxVarIndex = auxVarOffset + auxNum 

492 auxConIndex = auxConOffset + auxNum 

493 

494 # Prepare the auxiliary constraint. 

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

496 

497 if self._debug(): 

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

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

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

501 

502 # Add the auxiliary constraint. 

503 self.int.putarow(auxConIndex, J, V) 

504 self.int.putaij(auxConIndex, auxVarIndex, -1.0) 

505 for k, w in zip(K, W): 

506 self.int.putbaraij(auxConIndex, k, [w], [1.0]) 

507 self.int.putconbound(auxConIndex, mosek.boundkey.fx, -c, -c) 

508 

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

510 # variables (and auxiliary constraints). 

511 mosekVars[missingVarIndex] = auxVarIndex 

512 mosekCons[missingVarIndex] = auxConIndex 

513 

514 if self._debug(): 

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

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

517 

518 # Add the conic constraint. 

519 coneIndex = self.int.getnumcone() 

520 mosekCone = mosek.conetype.rquad if isRotated else mosek.conetype.quad 

521 self.int.appendcone(mosekCone, 0.0, mosekVars) 

522 

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

524 

525 def _import_sdp_constraint(self, picosConstraint): 

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

527 

528 import mosek 

529 

530 n = picosConstraint.size[0] 

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

532 

533 # MOSEK does not support general LMIs but so called "bar vars" which 

534 # are variables in the symmetric positive semidefinite cone. We use them 

535 # in combination with linear equalities to represent the LMI. 

536 barVar = self.int.getnumbarvar() 

537 mosekConOffset = self.int.getnumcon() 

538 self.int.appendbarvars([n]) 

539 self.int.appendcons(dim) 

540 

541 # MOSEK uses a storage of symmetric coefficient matrices that are used 

542 # as dot product coefficients to build scalar constraints involving both 

543 # "bar vars" and normal scalar variables. We build a couple of these 

544 # matrices to be able to select individual entries of our "bar vars". 

545 # More precisely, if X is a n×n symmetric matrix and 0 ≤ k ≤ n*(n+1)/2, 

546 # then <Units[n][k],X> = h(X)[k] where h refers to the lower-triangular, 

547 # column-major half-vectorization of X. 

548 if n in self._mosekBarUnitCoefs: 

549 Units = self._mosekBarUnitCoefs[n] 

550 else: 

551 Units = self._mosekBarUnitCoefs[n] = [ 

552 self.int.appendsparsesymmat( 

553 n, [row], [col], [1.0 if row == col else 0.5]) 

554 for col in range(n) for row in range(col, n)] 

555 

556 # We iterate over the lower triangular scalar sub-expressions of the 

557 # expression that the PICOS constraint states to be PSD, and constrain 

558 # them to be eqal to the MOSEK "bar var" at the same index. 

559 psdRows = tuple(self._affinexp_pic2msk(picosConstraint.psd)) 

560 for lowTriIndex, row, col in self._low_tri_indices(n): 

561 mosekConIndex = mosekConOffset + lowTriIndex 

562 J, V, K, W, c = psdRows[row + n*col] 

563 rhs = -c 

564 

565 # The lower-triangular entries in the PSD-constrained matrix … 

566 self.int.putarow(mosekConIndex, J, V) 

567 for k, w in zip(K, W): 

568 self.int.putbaraij(mosekConIndex, k, [w], [1.0]) 

569 

570 # … minus the corresponding bar var entries … 

571 self.int.putbaraij( 

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

573 

574 # … should equal zero. 

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

576 

577 if self._debug(): 

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

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

580 

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

582 

583 def _import_constraint(self, picosConstraint): 

584 if self._debug(): 

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

586 

587 if isinstance(picosConstraint, AffineConstraint): 

588 self._import_linear_constraint(picosConstraint) 

589 elif isinstance(picosConstraint, ConvexQuadraticConstraint): 

590 self._import_quad_constraint(picosConstraint) 

591 elif isinstance(picosConstraint, SOCConstraint) \ 

592 or isinstance(picosConstraint, RSOCConstraint): 

593 self._import_quad_conic_constraint(picosConstraint) 

594 elif isinstance(picosConstraint, LMIConstraint): 

595 self._import_sdp_constraint(picosConstraint) 

596 else: 

597 assert isinstance(picosConstraint, DummyConstraint), \ 

598 "Unexpected constraint type: {}".format( 

599 picosConstraint.__class__.__name__) 

600 

601 def _reset_objective(self): 

602 numVars = self.int.getnumvar() 

603 numBarVars = self.int.getnumbarvar() 

604 

605 # Reset affine part. 

606 self.int.putclist(range(numVars), [0.0]*numVars) 

607 for k in range(numBarVars): 

608 self.int.putbarcj(k, [], []) 

609 self.int.putcfix(0.0) 

610 

611 # Reset quadratic part. 

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

613 

614 def _import_affine_objective(self, picosObjective): 

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

616 

617 self.int.putclist(J, V) 

618 for k, w in zip(K, W): 

619 self.int.putbarcj(k, [w], [1.0]) 

620 self.int.putcfix(c) 

621 

622 def _import_quadratic_objective(self, picosObjective): 

623 # Mixing of conic and quadratic constraints is not supported by MOSEK; 

624 # therefor we do not handle bar vars here. 

625 assert not self._mosekBarVarIndex 

626 

627 # Import the quadratic part. 

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

629 

630 # Import the affine part. 

631 self._import_affine_objective(picosObjective.aff) 

632 

633 def _import_objective(self): 

634 import mosek 

635 

636 picosSense, picosObjective = self.ext.no 

637 

638 # Import objective sense. 

639 if picosSense == "min": 

640 self.int.putobjsense(mosek.objsense.minimize) 

641 else: 

642 assert picosSense == "max" 

643 self.int.putobjsense(mosek.objsense.maximize) 

644 

645 # Import objective function. 

646 if isinstance(picosObjective, AffineExpression): 

647 self._import_affine_objective(picosObjective) 

648 else: 

649 assert isinstance(picosObjective, QuadraticExpression) 

650 self._import_quadratic_objective(picosObjective) 

651 

652 def _import_problem(self): 

653 # Create a problem instance. 

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

655 

656 # Convert trivial LMIs and their symmetric variable into a PSD variable. 

657 # This allows the constraint to be loaded more efficiently. 

658 constraints, psdVars = [], set() 

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

660 if isinstance(constraint, LMIConstraint) and constraint.semidefVar: 

661 psdVars.add(constraint.semidefVar) 

662 else: 

663 constraints.append(constraint) 

664 

665 # Import normal variables. 

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

667 if variable not in psdVars: 

668 self._import_variable(variable) 

669 

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

671 for variable in psdVars: 

672 self._import_psd_variable(variable) 

673 

674 # Import remaining constraints. 

675 for constraint in constraints: 

676 self._import_constraint(constraint) 

677 

678 # Set objective. 

679 self._import_objective() 

680 

681 def _update_problem(self): 

682 for oldConstraint in self._removed_constraints(): 

683 raise ProblemUpdateError("PICOS does not support removing " 

684 "constraints from a MOSEK instance.") 

685 

686 for oldVariable in self._removed_variables(): 

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

688 " from a MOSEK instance.") 

689 

690 for newVariable in self._new_variables(): 

691 self._import_variable(newVariable) 

692 

693 for newConstraint in self._new_constraints(): 

694 self._import_constraint(newConstraint) 

695 

696 if self._objective_has_changed(): 

697 self._reset_objective() 

698 self._import_objective() 

699 

700 def _solve(self): 

701 import mosek 

702 

703 # Determine whether an LP is being solved. 

704 # TODO: Give Problem an interface for checks like this. 

705 _affine = (AffineConstraint, DummyConstraint) 

706 _is_lp = isinstance(self.ext.no.function, AffineExpression) \ 

707 and all(isinstance(constraint, _affine) 

708 for constraint in self.ext.constraints.values()) 

709 

710 # Reset all solver options to default. 

711 self.int.setdefaults() 

712 self.int.putoptserverhost("") 

713 

714 # verbosity 

715 if self._verbose(): 

716 self.int.set_Stream(mosek.streamtype.log, self._streamprinter) 

717 self.int.putintparam(mosek.iparam.log, self.ext.verbosity()) 

718 

719 # abs_prim_fsb_tol 

720 if self.ext.options.abs_prim_fsb_tol is not None: 

721 value = self.ext.options.abs_prim_fsb_tol 

722 

723 # Interior-point primal feasibility tolerances. 

724 self.int.putdouparam(mosek.dparam.intpnt_tol_pfeas, value) 

725 self.int.putdouparam(mosek.dparam.intpnt_co_tol_pfeas, value) 

726 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_pfeas, value) 

727 if self.ver <= 8: 

728 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_pfeas, value) 

729 

730 # Simplex primal feasibility tolerance. 

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

732 

733 # Mixed-integer (primal) feasibility tolerance. 

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

735 

736 # abs_dual_fsb_tol 

737 if self.ext.options.abs_dual_fsb_tol is not None: 

738 value = self.ext.options.abs_dual_fsb_tol 

739 

740 # Interior-point dual feasibility tolerances. 

741 self.int.putdouparam(mosek.dparam.intpnt_tol_dfeas, value) 

742 self.int.putdouparam(mosek.dparam.intpnt_co_tol_dfeas, value) 

743 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_dfeas, value) 

744 if self.ver <= 8: 

745 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_dfeas, value) 

746 

747 # Simplex dual feasibility (optimality) tolerance. 

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

749 

750 # rel_dual_fsb_tol 

751 if self.ext.options.rel_dual_fsb_tol is not None: 

752 # Simplex relative dual feasibility (optimality) tolerance. 

753 self.int.putdouparam(mosek.dparam.basis_rel_tol_s, 

754 self.ext.options.rel_dual_fsb_tol) 

755 

756 # rel_ipm_opt_tol 

757 if self.ext.options.rel_ipm_opt_tol is not None: 

758 value = self.ext.options.rel_ipm_opt_tol 

759 

760 # Interior-point primal feasibility tolerances. 

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

762 self.int.putdouparam(mosek.dparam.intpnt_co_tol_rel_gap, value) 

763 self.int.putdouparam(mosek.dparam.intpnt_qo_tol_rel_gap, value) 

764 if self.ver <= 8: 

765 self.int.putdouparam(mosek.dparam.intpnt_nl_tol_rel_gap, value) 

766 

767 # abs_bnb_opt_tol 

768 if self.ext.options.abs_bnb_opt_tol is not None: 

769 self.int.putdouparam(mosek.dparam.mio_tol_abs_gap, 

770 self.ext.options.abs_bnb_opt_tol) 

771 

772 # rel_bnb_opt_tol 

773 if self.ext.options.rel_bnb_opt_tol is not None: 

774 self.int.putdouparam(mosek.dparam.mio_tol_rel_gap, 

775 self.ext.options.rel_bnb_opt_tol) 

776 

777 # integrality_tol 

778 if self.ext.options.integrality_tol is not None: 

779 self.int.putdouparam(mosek.dparam.mio_tol_abs_relax_int, 

780 self.ext.options.integrality_tol) 

781 

782 # max_iterations 

783 if self.ext.options.max_iterations is not None: 

784 value = self.ext.options.max_iterations 

785 self.int.putintparam(mosek.iparam.bi_max_iterations, value) 

786 self.int.putintparam(mosek.iparam.intpnt_max_iterations, value) 

787 self.int.putintparam(mosek.iparam.sim_max_iterations, value) 

788 

789 # Prepare lp_node_method and lp_root_method. 

790 _lpm = { 

791 "interior": (mosek.optimizertype.intpnt if _is_lp 

792 else mosek.optimizertype.conic), 

793 "psimplex": mosek.optimizertype.primal_simplex, 

794 "dsimplex": mosek.optimizertype.dual_simplex} 

795 

796 # lp_node_method 

797 if self.ext.options.lp_node_method is not None: 

798 value = self.ext.options.lp_node_method 

799 assert value in _lpm, "Unexpected lp_node_method value." 

800 self.int.putintparam(mosek.iparam.mio_node_optimizer, _lpm[value]) 

801 

802 # lp_root_method 

803 if self.ext.options.lp_root_method is not None: 

804 value = self.ext.options.lp_root_method 

805 assert value in _lpm, "Unexpected lp_root_method value." 

806 self.int.putintparam(mosek.iparam.mio_root_optimizer, _lpm[value]) 

807 

808 # timelimit 

809 if self.ext.options.timelimit is not None: 

810 value = float(self.ext.options.timelimit) 

811 self.int.putdouparam(mosek.dparam.optimizer_max_time, value) 

812 self.int.putdouparam(mosek.dparam.mio_max_time, value) 

813 

814 # max_fsb_nodes 

815 if self.ext.options.max_fsb_nodes is not None: 

816 self.int.putintparam(mosek.iparam.mio_max_num_solutions, 

817 self.ext.options.max_fsb_nodes) 

818 

819 # Handle MOSEK-specific options. 

820 for key, value in self.ext.options.mosek_params.items(): 

821 try: 

822 self.int.putparam(key.upper(), str(value)) 

823 except mosek.Error as error: 

824 self._handle_bad_solver_specific_option(key, value, error) 

825 

826 # Handle 'mosek_basic_sol' option. 

827 if self.ext.options.mosek_basic_sol and _is_lp: 

828 _intpnt_basis = mosek.basindtype.always 

829 else: 

830 _intpnt_basis = mosek.basindtype.never 

831 self.int.putintparam(mosek.iparam.intpnt_basis, _intpnt_basis) 

832 

833 # Handle 'mosek_server' option. 

834 if self.ext.options.mosek_server: 

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

836 

837 # Handle unsupported options. 

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

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

840 

841 # Attempt to solve the problem. 

842 with self._header(), self._stopwatch(): 

843 try: 

844 self.int.optimize() 

845 except mosek.Error as error: 

846 if error.errno in ( 

847 mosek.rescode.err_con_q_not_psd, 

848 mosek.rescode.err_con_q_not_nsd, 

849 mosek.rescode.err_obj_q_not_psd, 

850 mosek.rescode.err_obj_q_not_nsd, 

851 mosek.rescode.err_toconic_constr_q_not_psd, 

852 mosek.rescode.err_toconic_objective_not_psd): 

853 self._handle_continuous_nonconvex_error(error) 

854 else: 

855 raise 

856 

857 # Set the solution to be retrieved. 

858 if self.ext.is_continuous(): 

859 if _intpnt_basis == mosek.basindtype.always: 

860 solType = mosek.soltype.bas 

861 else: 

862 solType = mosek.soltype.itr 

863 else: 

864 solType = mosek.soltype.itg 

865 

866 # Retrieve primals. 

867 primals = {} 

868 if self.ext.options.primals is not False: 

869 values = [float("nan")]*self.int.getnumvar() 

870 self.int.getxx(solType, values) 

871 

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

873 if picosVar in self._mosekBarVarIndex: 

874 barVarIndex = self._mosekBarVarIndex[picosVar] 

875 

876 assert isinstance(picosVar, SymmetricVariable) 

877 

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

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

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

881 

882 # Obtain a lower triangular vectorization from MOSEK. 

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

884 

885 # Convert to a full vectorization. 

886 # TODO: Immediately convert to a symmetric vectorization. 

887 for lti, row, col in self._low_tri_indices(n): 

888 value = lowTriPrimal[lti] 

889 matrixPrimal[n*row + col] = value 

890 if row != col: 

891 matrixPrimal[n*col + row] = value 

892 

893 # Load the full vectorization as a CVXOPT matrix. 

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

895 

896 # Re-vectorize using a symmetric vectorization. 

897 primal = picosVar._vec.vectorize(matrixPrimal) 

898 else: 

899 mosekOffset = self._mosekVarOffset[picosVar] 

900 primal = values[mosekOffset:mosekOffset + picosVar.dim] 

901 

902 if float("nan") in primal: 

903 primals[picosVar] = None 

904 else: 

905 primals[picosVar] = primal 

906 

907 # Retrieve duals. 

908 duals = {} 

909 if self.ext.options.duals is not False and self.ext.is_continuous(): 

910 minimizing = self.int.getobjsense() == mosek.objsense.minimize 

911 

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

913 if isinstance(constraint, DummyConstraint): 

914 duals[constraint] = cvxopt.spmatrix( 

915 [], [], [], constraint.size) 

916 continue 

917 

918 length = len(constraint) 

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

920 

921 if isinstance(constraint, AffineConstraint): 

922 offset = self._mosekLinConOffset[constraint] 

923 self.int.getyslice(solType, offset, offset + length, dual) 

924 elif isinstance(constraint, ConvexQuadraticConstraint): 

925 # TODO: Implement consistent QCQP dual retrieval for all 

926 # solvers that return such duals. 

927 dual = None 

928 elif isinstance(constraint, SOCConstraint) \ 

929 or isinstance(constraint, RSOCConstraint): 

930 mosekVars = self._mosekCone[constraint][1] 

931 for localConeIndex in range(length): 

932 x = [float("nan")] 

933 offset = mosekVars[localConeIndex] 

934 self.int.getsnxslice(solType, offset, offset + 1, x) 

935 dual[localConeIndex] = x[0] 

936 

937 if isinstance(constraint, SOCConstraint): 

938 dual = [-du for du in dual] 

939 elif isinstance(constraint, RSOCConstraint): 

940 dual[0] = -dual[0] / math.sqrt(2.0) 

941 dual[1] = -dual[1] / math.sqrt(2.0) 

942 dual[2:] = [-du for du in dual[2:]] 

943 elif isinstance(constraint, LMIConstraint): 

944 if constraint.semidefVar: 

945 assert constraint not in self._mosekLMI 

946 picosVar = constraint.semidefVar 

947 assert picosVar in self._mosekBarVarIndex 

948 barVar = self._mosekBarVarIndex[picosVar] 

949 else: 

950 assert not constraint.semidefVar 

951 assert constraint in self._mosekLMI 

952 barVar, _ = self._mosekLMI[constraint] 

953 

954 n = constraint.size[0] 

955 d = n*(n + 1) // 2 

956 lowerTriangularDual = [float("nan")]*d 

957 self.int.getbarsj(solType, barVar, lowerTriangularDual) 

958 for lti, row, col in self._low_tri_indices(n): 

959 value = lowerTriangularDual[lti] 

960 dual[n*row + col] = value 

961 if row != col: 

962 dual[n*col + row] = value 

963 else: 

964 assert False, "Constraint type not supported." 

965 

966 if dual is None: 

967 pass 

968 elif float("nan") in dual: 

969 dual = None 

970 else: 

971 dual = cvxopt.matrix(dual, constraint.size) 

972 

973 if isinstance(constraint, AffineConstraint) \ 

974 or isinstance(constraint, LMIConstraint): 

975 if not constraint.is_increasing(): 

976 dual = -dual 

977 

978 if minimizing: 

979 dual = -dual 

980 

981 duals[constraint] = dual 

982 

983 # Retrieve objective value. 

984 value = self.int.getprimalobj(solType) 

985 

986 # Retrieve solution status. 

987 primalStatus = self._solution_status_pic2msk( 

988 self.int.getsolsta(solType), primalSolution=True) 

989 dualStatus = self._solution_status_pic2msk( 

990 self.int.getsolsta(solType), primalSolution=False) 

991 problemStatus = self._problem_status_pic2msk( 

992 self.int.getprosta(solType), not self.ext.is_continuous()) 

993 

994 return self._make_solution( 

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

996 

997 def _solution_status_pic2msk(self, statusCode, primalSolution): 

998 from mosek import solsta as ss 

999 dualSolution = not primalSolution 

1000 

1001 map = { 

1002 ss.unknown: SS_UNKNOWN, 

1003 ss.optimal: SS_OPTIMAL, 

1004 ss.prim_feas: 

1005 SS_FEASIBLE if primalSolution else SS_UNKNOWN, 

1006 ss.dual_feas: 

1007 SS_FEASIBLE if dualSolution else SS_UNKNOWN, 

1008 ss.prim_and_dual_feas: SS_FEASIBLE, 

1009 ss.prim_infeas_cer: 

1010 SS_INFEASIBLE if primalSolution else SS_UNKNOWN, 

1011 ss.dual_infeas_cer: 

1012 SS_INFEASIBLE if dualSolution else SS_UNKNOWN, 

1013 ss.prim_illposed_cer: SS_UNKNOWN, 

1014 ss.dual_illposed_cer: SS_UNKNOWN, 

1015 ss.integer_optimal: SS_OPTIMAL, 

1016 } 

1017 

1018 if self.ver < 9: 

1019 map.update({ 

1020 ss.near_optimal: SS_FEASIBLE, 

1021 ss.near_prim_feas: SS_UNKNOWN, 

1022 ss.near_dual_feas: SS_UNKNOWN, 

1023 ss.near_prim_and_dual_feas: SS_UNKNOWN, 

1024 ss.near_prim_infeas_cer: SS_UNKNOWN, 

1025 ss.near_dual_infeas_cer: SS_UNKNOWN, 

1026 ss.near_integer_optimal: SS_FEASIBLE 

1027 }) 

1028 

1029 try: 

1030 return map[statusCode] 

1031 except KeyError: 

1032 self._warn( 

1033 "The MOSEK solution status code {} is not known to PICOS." 

1034 .format(statusCode)) 

1035 return SS_UNKNOWN 

1036 

1037 def _problem_status_pic2msk(self, statusCode, integerProblem): 

1038 from mosek import prosta as ps 

1039 

1040 map = { 

1041 ps.unknown: PS_UNKNOWN, 

1042 ps.prim_and_dual_feas: PS_FEASIBLE, 

1043 ps.prim_feas: 

1044 PS_FEASIBLE if integerProblem else PS_UNKNOWN, 

1045 ps.dual_feas: PS_UNKNOWN, 

1046 ps.prim_infeas: PS_INFEASIBLE, 

1047 ps.dual_infeas: PS_UNBOUNDED, # TODO: UNB_OR_INF 

1048 ps.prim_and_dual_infeas: PS_INFEASIBLE, 

1049 ps.ill_posed: PS_ILLPOSED, 

1050 ps.prim_infeas_or_unbounded: PS_UNKNOWN # TODO: UNB_OR_INF 

1051 } 

1052 

1053 if self.ver < 9: 

1054 map.update({ 

1055 ps.near_prim_and_dual_feas: PS_UNKNOWN, 

1056 ps.near_prim_feas: PS_UNKNOWN, 

1057 ps.near_dual_feas: PS_UNKNOWN 

1058 }) 

1059 

1060 try: 

1061 return map[statusCode] 

1062 except KeyError: 

1063 self._warn("The MOSEK problem status code {} is not known to PICOS." 

1064 .format(statusCode)) 

1065 return PS_UNKNOWN 

1066 

1067 

1068# -------------------------------------- 

1069__all__ = api_end(_API_START, globals())