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) 

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._mosekLinConOffset = dict() 

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

139 constraint (start) index.""" 

140 

141 self._mosekQuadConIndex = dict() 

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

143 

144 self._mosekCone = dict() 

145 """ 

146 Maps a PICOS SOCC or RSOCC to a triple: 

147 

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

159 

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

164 

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

172 

173 def reset_problem(self): 

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

175 self.int = None 

176 

177 self._mosekVarOffset.clear() 

178 self._mosekLinConOffset.clear() 

179 self._mosekQuadConIndex.clear() 

180 self._mosekCone.clear() 

181 self._mosekLMI.clear() 

182 self._mosekBarUnitCoefs.clear() 

183 

184 @classmethod 

185 def _get_environment(cls): 

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

187 import mosek 

188 cls.mosekEnvironment = mosek.Env() 

189 

190 return cls.mosekEnvironment 

191 

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

196 

197 @classmethod 

198 def _get_major_version(cls): 

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

200 import mosek 

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

202 

203 return cls.mosekVersion[0] 

204 

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

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

207 

208 @staticmethod 

209 def _streamprinter(text): 

210 sys.stdout.write(text) 

211 sys.stdout.flush() 

212 

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) 

219 

220 def _scalar_affinexp_pic2msk(self, picosExpression): 

221 assert len(picosExpression) == 1 

222 return picosExpression.sparse_rows(self._mosekVarOffset)[0] 

223 

224 def _affinexp_pic2msk(self, picosExpression): 

225 return picosExpression.sparse_rows(self._mosekVarOffset) 

226 

227 def _quadexp_pic2msk(self, picosExpression): 

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

229 

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

238 

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 

252 

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 

259 

260 I.append(mskVar1Index) 

261 J.append(mskVar2Index) 

262 V.append(localCoefficient) 

263 

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 

269 

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

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

272 

273 def _import_variable(self, picosVar): 

274 import mosek 

275 

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

277 dim = picosVar.dim 

278 indices = range(numVars, numVars + dim) 

279 self.int.appendvars(dim) 

280 

281 # Set the variable type. 

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

283 self.int.putvartypelist( 

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

285 

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 

294 

295 lower, upper = picosVar.bound_dicts 

296 

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

298 boundKeys[i] = mosek.boundkey.lo 

299 lowerBounds[i] = b 

300 

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 

309 

310 upperBounds[i] = b 

311 

312 self.int.putvarboundlist( 

313 indices, boundKeys, lowerBounds, upperBounds) 

314 

315 def _import_linear_constraint(self, picosConstraint): 

316 import mosek 

317 

318 numCons = self.int.getnumcon() 

319 conLen = len(picosConstraint) 

320 self.int.appendcons(conLen) 

321 

322 rows = picosConstraint.sparse_Ab_rows(self._mosekVarOffset) 

323 

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 

330 

331 for localConIndex, (mosekVarIndices, coefs, offset) in enumerate(rows): 

332 mosekConIndex = numCons + localConIndex 

333 

334 self.int.putarow(mosekConIndex, mosekVarIndices, coefs) 

335 self.int.putconbound(mosekConIndex, boundKey, offset, offset) 

336 

337 self._mosekLinConOffset[picosConstraint] = numCons 

338 

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) 

344 

345 # Add the quadratic part. 

346 self.int.putqconk( 

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

348 

349 self._mosekQuadConIndex[picosConstraint] = mosekConIndex 

350 

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 

358 

359 def _import_quad_conic_constraint(self, picosConstraint): 

360 import mosek 

361 

362 isRotated = isinstance(picosConstraint, RSOCConstraint) 

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

364 

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) 

377 

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) 

387 

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) 

394 

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) 

398 

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 

404 

405 # Prepare the auxiliary constraint. 

406 indices, values, offset = entries[missingVarIndex] 

407 

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) 

420 

421 if self._debug(): 

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

423 "{}.T * x{} = {}".format(values, indices, -offset)) 

424 

425 # Add the auxiliary constraint. 

426 self.int.putarow(auxConIndex, indices, values) 

427 self.int.putconbound( 

428 auxConIndex, mosek.boundkey.fx, -offset, -offset) 

429 

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

431 # variables (and auxiliary constraints). 

432 mosekVars[missingVarIndex] = auxVarIndex 

433 mosekCons[missingVarIndex] = auxConIndex 

434 

435 if self._debug(): 

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

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

438 

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) 

443 

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

445 

446 def _import_sdp_constraint(self, picosConstraint): 

447 import mosek 

448 

449 n = picosConstraint.size[0] 

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

451 

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) 

459 

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

474 

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] 

483 

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

485 self.int.putarow(mosekConIndex, indices, coefficients) 

486 

487 # … minus the corresponding bar var entries … 

488 self.int.putbaraij( 

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

490 

491 # … should equal zero. 

492 self.int.putconbound( 

493 mosekConIndex, mosek.boundkey.fx, -offset, -offset) 

494 

495 if self._debug(): 

496 self._debug(" Index {} ({}, {}): indices = {}, coefs = {}" 

497 .format(lowTriIndex, row, col, indices, coefficients)) 

498 

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

500 

501 def _import_constraint(self, picosConstraint): 

502 if self._debug(): 

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

504 

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

518 

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) 

524 

525 # Reset quadratic part. 

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

527 

528 def _import_affine_objective(self, picosObjective): 

529 mosekIndices = [] 

530 mosekCoefficients = [] 

531 

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

538 

539 self.int.putclist(mosekIndices, mosekCoefficients) 

540 

541 if picosObjective._constant_coef: 

542 self.int.putcfix(picosObjective._constant_coef[0]) 

543 

544 def _import_quadratic_objective(self, picosObjective): 

545 # Import the quadratic part. 

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

547 

548 # Import the affine part. 

549 self._import_affine_objective(picosObjective.aff) 

550 

551 def _import_objective(self): 

552 import mosek 

553 

554 picosSense, picosObjective = self.ext.no 

555 

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) 

562 

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) 

569 

570 def _import_problem(self): 

571 # Create a problem instance. 

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

573 

574 # Import variables. 

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

576 self._import_variable(variable) 

577 

578 # Import constraints. 

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

580 self._import_constraint(constraint) 

581 

582 # Set objective. 

583 self._import_objective() 

584 

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

589 

590 for oldVariable in self._removed_variables(): 

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

592 " from a MOSEK instance.") 

593 

594 for newVariable in self._new_variables(): 

595 self._import_variable(newVariable) 

596 

597 for newConstraint in self._new_constraints(): 

598 self._import_constraint(newConstraint) 

599 

600 if self._objective_has_changed(): 

601 self._reset_objective() 

602 self._import_objective() 

603 

604 def _solve(self): 

605 import mosek 

606 

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

612 

613 # Reset all solver options to default. 

614 self.int.setdefaults() 

615 

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

620 

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 

624 

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) 

631 

632 # Simplex primal feasibility tolerance. 

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

634 

635 # Mixed-integer (primal) feasibility tolerance. 

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

637 

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 

641 

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) 

648 

649 # Simplex dual feasibility (optimality) tolerance. 

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

651 

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) 

657 

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 

661 

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) 

668 

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) 

673 

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) 

678 

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) 

683 

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) 

690 

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} 

696 

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

702 

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

708 

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) 

714 

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) 

719 

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) 

726 

727 # Handle unsupported options. 

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

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

730 

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 

746 

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 

755 

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) 

761 

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

763 mosekOffset = self._mosekVarOffset[picosVar] 

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

765 

766 if float("nan") in primal: 

767 primals[picosVar] = None 

768 else: 

769 primals[picosVar] = primal 

770 

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 

775 

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

777 if isinstance(constraint, DummyConstraint): 

778 duals[constraint] = cvxopt.spmatrix( 

779 [], [], [], constraint.size) 

780 continue 

781 

782 length = len(constraint) 

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

784 

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] 

800 

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

819 

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) 

826 

827 if isinstance(constraint, AffineConstraint) \ 

828 or isinstance(constraint, LMIConstraint): 

829 if not constraint.is_increasing(): 

830 dual = -dual 

831 

832 if minimizing: 

833 dual = -dual 

834 

835 duals[constraint] = dual 

836 

837 # Retrieve objective value. 

838 value = self.int.getprimalobj(solType) 

839 

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

847 

848 return self._make_solution( 

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

850 

851 def _solution_status_pic2msk(self, statusCode, primalSolution): 

852 from mosek import solsta as ss 

853 dualSolution = not primalSolution 

854 

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 } 

871 

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

882 

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 

890 

891 def _problem_status_pic2msk(self, statusCode, integerProblem): 

892 from mosek import prosta as ps 

893 

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 } 

906 

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

913 

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 

920 

921 

922# -------------------------------------- 

923__all__ = api_end(_API_START, globals())