Coverage for picos/solvers/solver_gurobi.py: 77.55%

548 statements  

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

1# ------------------------------------------------------------------------------ 

2# Copyright (C) 2018-2022 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:`GurobiSolver`.""" 

20 

21from collections import namedtuple 

22 

23import cvxopt 

24import numpy 

25 

26from .. import settings 

27from ..apidoc import api_end, api_start 

28from ..constraints import (AffineConstraint, ConvexQuadraticConstraint, 

29 DummyConstraint, NonconvexQuadraticConstraint, 

30 RSOCConstraint, SOCConstraint) 

31from ..expressions import (CONTINUOUS_VARTYPES, AffineExpression, 

32 BinaryVariable, IntegerVariable, 

33 QuadraticExpression) 

34from ..expressions.data import cvx2csr, cvx2np 

35from ..modeling.footprint import Specification 

36from ..modeling.solution import (PS_FEASIBLE, PS_INF_OR_UNB, PS_INFEASIBLE, 

37 PS_UNBOUNDED, PS_UNKNOWN, PS_UNSTABLE, 

38 SS_EMPTY, SS_FEASIBLE, SS_INFEASIBLE, 

39 SS_OPTIMAL, SS_PREMATURE, SS_UNKNOWN) 

40from .solver import Solver 

41 

42_API_START = api_start(globals()) 

43# ------------------------------- 

44 

45 

46class GurobiSolver(Solver): 

47 """Interface to the Gurobi solver via its official Python interface.""" 

48 

49 # TODO: Don't support (conic) quadratic constraints when duals are 

50 # requested because their precision is bad and can't be controlled? 

51 SUPPORTED_8 = Specification( 

52 objectives=[ 

53 AffineExpression, 

54 QuadraticExpression], 

55 constraints=[ 

56 DummyConstraint, 

57 AffineConstraint, 

58 SOCConstraint, 

59 RSOCConstraint, 

60 ConvexQuadraticConstraint]) 

61 

62 SUPPORTED_9 = Specification( 

63 objectives=[ 

64 AffineExpression, 

65 QuadraticExpression], 

66 constraints=[ 

67 DummyConstraint, 

68 AffineConstraint, 

69 SOCConstraint, 

70 RSOCConstraint, 

71 NonconvexQuadraticConstraint]) 

72 

73 @classmethod 

74 def _gurobi9(cls): 

75 try: 

76 import gurobipy as gurobi 

77 except ImportError: 

78 # This method should be used only after test_availability confirmed 

79 # that gurobipy is available, however that method does not actually 

80 # perform an import and so it could still fail here due to a bad 

81 # installation. In this case an exception will be raised when Gurobi 

82 # is actually selected and we just return False as a dummy here. 

83 return False 

84 else: 

85 return hasattr(gurobi, "MVar") 

86 

87 @classmethod 

88 def supports(cls, footprint, explain=False): 

89 """Implement :meth:`~.solver.Solver.supports`.""" 

90 result = Solver.supports(footprint, explain) 

91 if not result or (explain and not result[0]): 

92 return result 

93 

94 supported = cls.SUPPORTED_9 if cls._gurobi9() else cls.SUPPORTED_8 

95 

96 if footprint not in supported: 

97 if explain: 

98 return False, supported.mismatch_reason(footprint) 

99 else: 

100 return False 

101 

102 return (True, None) if explain else True 

103 

104 @classmethod 

105 def default_penalty(cls): 

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

107 return 0.0 # Commercial solver. 

108 

109 @classmethod 

110 def test_availability(cls): 

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

112 cls.check_import("gurobipy") 

113 

114 @classmethod 

115 def names(cls): 

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

117 return "gurobi", "Gurobi", "Gurobi Optimizer", None 

118 

119 @classmethod 

120 def is_free(cls): 

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

122 return False 

123 

124 GurobiMetaConstraint = namedtuple( 

125 "GurobiMetaConstraint", ("auxCons", "auxVars")) 

126 

127 def __init__(self, problem): 

128 """Initialize a Gurobi solver interface. 

129 

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

131 """ 

132 super(GurobiSolver, self).__init__(problem) 

133 

134 self._matint_decision = None 

135 

136 self._gurobiVar = dict() 

137 """Maps PICOS variable indices to Gurobi variables. (Matrix interf.)""" 

138 

139 self._gurobiVars = [] 

140 """A list of all Gurobi variables added. (Legacy interface.)""" 

141 

142 self._gurobiVarOffset = dict() 

143 """Maps PICOS variables to Gurobi variable offsets. (Legacy interf.)""" 

144 

145 self._gurobiLinCon = dict() 

146 """Maps a PICOS linear constraint to (a) Gurobi linear constraint(s).""" 

147 

148 self._gurobiQuadCon = dict() 

149 """Maps a PICOS quadr. constraint to (a) Gurobi quadr. constraint(s).""" 

150 

151 self._gurobiConicCon = dict() 

152 """Maps a PICOS quadr. constraint to its Gurobi representation.""" 

153 

154 def _make_matint_decision(self): 

155 default = settings.PREFER_GUROBI_MATRIX_INTERFACE 

156 choice = self.ext.options.gurobi_matint 

157 

158 if not choice and (not default or choice is not None): 

159 return False 

160 

161 if not self._gurobi9(): 

162 if choice: 

163 raise RuntimeError( 

164 "Gurobi's matrix interface should be used by user's choice " 

165 "but Gurobi < 9 appears to be installed. More precisely, " 

166 "gurobipy.Mvar is not available.") 

167 else: 

168 return False 

169 

170 try: 

171 self.check_import("scipy.sparse") 

172 except ModuleNotFoundError as error: 

173 if choice: 

174 raise RuntimeError( 

175 "Gurobi's matrix interface should be used by user's choice " 

176 "but this requires SciPy, which was not found.") from error 

177 else: 

178 return False 

179 

180 assert choice or (choice is None and default) 

181 return True 

182 

183 @property 

184 def matint(self): 

185 """Whether Gurobi's matrix interface is in use.""" 

186 if self._matint_decision is None: 

187 self._matint_decision = self._make_matint_decision() 

188 

189 decision = self._matint_decision 

190 choice = self.ext.options.gurobi_matint 

191 

192 if (choice and not decision) \ 

193 or (decision and not choice and choice is not None): 

194 raise NotImplementedError( 

195 "The user's choice with respect to using Gurobi's matrix " 

196 "interface has changed between solution attempts. This is not " 

197 "supported. To re-load the problem with the other interface, " 

198 "you must manually reset your problem's solution strategy.") 

199 

200 return decision 

201 

202 def reset_problem(self): 

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

204 self.int = None 

205 

206 self._gurobiVar.clear() 

207 self._gurobiVars.clear() 

208 self._gurobiVarOffset.clear() 

209 

210 self._gurobiLinCon.clear() 

211 self._gurobiQuadCon.clear() 

212 self._gurobiConicCon.clear() 

213 

214 def _import_variable(self, picosVar): 

215 import gurobipy as gurobi 

216 

217 dim = picosVar.dim 

218 

219 # Retrieve types. 

220 if isinstance(picosVar, CONTINUOUS_VARTYPES): 

221 gurobiVarType = gurobi.GRB.CONTINUOUS 

222 elif isinstance(picosVar, IntegerVariable): 

223 gurobiVarType = gurobi.GRB.INTEGER 

224 elif isinstance(picosVar, BinaryVariable): 

225 gurobiVarType = gurobi.GRB.BINARY 

226 else: 

227 assert False, "Unexpected variable type." 

228 

229 # Retrieve bounds. 

230 lowerBounds = [-gurobi.GRB.INFINITY]*dim 

231 upperBounds = [gurobi.GRB.INFINITY]*dim 

232 lower, upper = picosVar.bound_dicts 

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

234 lowerBounds[i] = b 

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

236 upperBounds[i] = b 

237 

238 # Import the variable. 

239 if self.matint: 

240 gurobiVar = self.int.addMVar(dim, lb=lowerBounds, ub=upperBounds, 

241 vtype=gurobiVarType, name=picosVar.name) 

242 

243 self._gurobiVar[picosVar] = gurobiVar 

244 else: 

245 gurobiVarsDict = self.int.addVars( 

246 dim, lb=lowerBounds, ub=upperBounds, vtype=gurobiVarType) 

247 gurobiVars = [gurobiVarsDict[i] for i in range(dim)] 

248 

249 self._gurobiVarOffset[picosVar] = len(self._gurobiVars) 

250 self._gurobiVars.extend(gurobiVars) 

251 

252 def _remove_variable(self, picosVar): 

253 if self.matint: 

254 gurobiVar = self._gurobiVar.pop(picosVar) 

255 

256 self.int.remove(gurobiVar) 

257 else: 

258 offset = self._gurobiVarOffset[picosVar] 

259 dim = picosVar.dim 

260 

261 gurobiVars = self._gurobiVars[offset:offset + dim] 

262 

263 self._gurobiVars = ( 

264 self._gurobiVars[:offset] + self._gurobiVars[offset + dim:]) 

265 

266 for other in self._gurobiVarOffset: 

267 if self._gurobiVarOffset[other] > offset: 

268 self._gurobiVarOffset[other] -= dim 

269 

270 self.int.remove(gurobiVars) 

271 

272 def _import_variable_values(self): 

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

274 if picosVar.valued: 

275 value = picosVar.internal_value 

276 

277 if self.matint: 

278 gurobiVar = self._gurobiVar[picosVar] 

279 gurobiVar.Start = value 

280 else: 

281 offset = self._gurobiVarOffset[picosVar] 

282 dim = picosVar.dim 

283 

284 gurobiVars = self._gurobiVars[offset, offset + picosVar.dim] 

285 

286 for localIndex in range(dim): 

287 gurobiVars[localIndex].Start = value[localIndex] 

288 

289 def _reset_variable_values(self): 

290 import gurobipy as gurobi 

291 

292 if self.matint: 

293 gurobiVars = self._gurobiVar.values() 

294 else: 

295 gurobiVars = self._gurobiVars 

296 

297 for gurobiVar in gurobiVars: 

298 gurobiVar.Start = gurobi.GRB.UNDEFINED 

299 

300 def _affexp_pic2grb_matint(self, picosExpression): 

301 assert self.matint 

302 

303 # NOTE: Constant Gurobi matrix expressions don't exist; return thus 

304 # constant PICOS expressions as NumPy arrays. 

305 gurobiExpression = numpy.ravel(cvx2np( 

306 picosExpression._constant_coef)) 

307 

308 for picosVar, coef in picosExpression._linear_coefs.items(): 

309 A = cvx2csr(coef) 

310 x = self._gurobiVar[picosVar] 

311 

312 # NOTE: Using __(r)matmul__ as PICOS supports Python 3.4 and the 

313 # @-operator was implemented in Python 3.5. 

314 gurobiExpression += x.__rmatmul__(A) 

315 

316 return gurobiExpression 

317 

318 def _affexp_pic2grb_legacy(self, picosExpression): 

319 import gurobipy as gurobi 

320 

321 assert not self.matint 

322 

323 for J, V, c in picosExpression.sparse_rows(self._gurobiVarOffset): 

324 gurobiVars = [self._gurobiVars[j] for j in J] 

325 gurobiExpression = gurobi.LinExpr(V, gurobiVars) 

326 gurobiExpression.addConstant(c) 

327 

328 yield gurobiExpression 

329 

330 def _scalar_affexp_pic2grb(self, picosExpression): 

331 assert len(picosExpression) == 1 

332 

333 if self.matint: 

334 gurobiExpression = self._affexp_pic2grb_matint(picosExpression) 

335 

336 if picosExpression.constant: 

337 assert isinstance(gurobiExpression, numpy.ndarray) 

338 assert gurobiExpression.shape == (1,) 

339 

340 return gurobiExpression[0] 

341 else: 

342 return gurobiExpression 

343 else: 

344 return next(self._affexp_pic2grb_legacy(picosExpression)) 

345 

346 def _quadexp_pic2grb(self, picosExpression): 

347 import gurobipy as gurobi 

348 

349 assert isinstance(picosExpression, QuadraticExpression) 

350 

351 if self.matint: 

352 # Import affine part. 

353 gurobiExpression = self._affexp_pic2grb_matint(picosExpression.aff) 

354 

355 # Import quadratic part. 

356 for picosVars, coef in picosExpression._sparse_quads.items(): 

357 Q = cvx2csr(coef) 

358 x = self._gurobiVar[picosVars[0]] 

359 y = self._gurobiVar[picosVars[1]] 

360 

361 gurobiExpression += x.__matmul__(Q).__matmul__(y) 

362 else: 

363 # Import affine part. 

364 gurobiExpression = gurobi.QuadExpr( 

365 self._scalar_affexp_pic2grb(picosExpression.aff)) 

366 

367 # Import quadratic part. 

368 V, I, J = [], [], [] 

369 for (x, y), Q in picosExpression._sparse_quads.items(): 

370 xOffset = self._gurobiVarOffset[x] 

371 yOffset = self._gurobiVarOffset[y] 

372 

373 V.extend(Q.V) 

374 I.extend(self._gurobiVars[i] for i in Q.I + xOffset) 

375 J.extend(self._gurobiVars[j] for j in Q.J + yOffset) 

376 

377 gurobiExpression.addTerms(V, I, J) 

378 

379 return gurobiExpression 

380 

381 def _import_linear_constraint(self, picosCon): 

382 import gurobipy as gurobi 

383 

384 assert isinstance(picosCon, AffineConstraint) 

385 

386 if self.matint: 

387 gurobiLHS = self._affexp_pic2grb_matint(picosCon.lhs) 

388 gurobiRHS = self._affexp_pic2grb_matint(picosCon.rhs) 

389 

390 # HACK: Fallback to the legacy interface for constant constraints. 

391 # NOTE: This happens to work with remove_constraint since 

392 # gurobipy.Model.remove accepts both lists and constraints. 

393 if isinstance(gurobiLHS, numpy.ndarray) \ 

394 and isinstance(gurobiRHS, numpy.ndarray): 

395 if picosCon.is_increasing(): 

396 gurobiSense = gurobi.GRB.LESS_EQUAL 

397 elif picosCon.is_decreasing(): 

398 gurobiSense = gurobi.GRB.GREATER_EQUAL 

399 elif picosCon.is_equality(): 

400 gurobiSense = gurobi.GRB.EQUAL 

401 else: 

402 assert False, "Unexpected constraint relation." 

403 

404 return [self.int.addLConstr(a, gurobiSense, b) 

405 for a, b in zip(gurobiLHS, gurobiRHS)] 

406 

407 # Construct the constraint. 

408 if picosCon.is_increasing(): 

409 gurobiCon = gurobiLHS <= gurobiRHS 

410 elif picosCon.is_decreasing(): 

411 gurobiCon = gurobiLHS >= gurobiRHS 

412 elif picosCon.is_equality(): 

413 gurobiCon = gurobiLHS == gurobiRHS 

414 else: 

415 assert False, "Unexpected constraint relation." 

416 

417 # Add the constraint. 

418 gurobiCon = self.int.addConstr(gurobiCon) 

419 

420 return gurobiCon 

421 else: 

422 # Retrieve sense. 

423 if picosCon.is_increasing(): 

424 gurobiSense = gurobi.GRB.LESS_EQUAL 

425 elif picosCon.is_decreasing(): 

426 gurobiSense = gurobi.GRB.GREATER_EQUAL 

427 elif picosCon.is_equality(): 

428 gurobiSense = gurobi.GRB.EQUAL 

429 else: 

430 assert False, "Unexpected constraint relation." 

431 

432 # Append scalar constraints. 

433 gurobiCons = [self.int.addLConstr(gurobiLHS, gurobiSense, 0.0) 

434 for gurobiLHS in self._affexp_pic2grb_legacy(picosCon.lmr)] 

435 

436 return gurobiCons 

437 

438 def _import_quad_constraint(self, picosCon): 

439 import gurobipy as gurobi 

440 

441 # NOTE: NonconvexQuadraticConstraint includes ConvexQuadraticConstraint. 

442 assert isinstance(picosCon, NonconvexQuadraticConstraint) 

443 

444 if self.matint: 

445 gurobiLE0 = self._quadexp_pic2grb(picosCon.le0) 

446 gurobiCon = self.int.addConstr(gurobiLE0 <= 0) 

447 else: 

448 gurobiLHS = self._quadexp_pic2grb(picosCon.le0) 

449 gurobiRHS = -gurobiLHS.getLinExpr().getConstant() 

450 

451 if gurobiRHS: 

452 gurobiLHS.getLinExpr().addConstant(gurobiRHS) 

453 

454 gurobiCon = self.int.addQConstr( 

455 gurobiLHS, gurobi.GRB.LESS_EQUAL, gurobiRHS) 

456 

457 return gurobiCon 

458 

459 def _import_socone_constraint(self, picosCon): 

460 import gurobipy as gurobi 

461 

462 assert isinstance(picosCon, SOCConstraint) 

463 

464 n = len(picosCon.ne) 

465 

466 # Load defining expressions. 

467 gurobiRHS = self._scalar_affexp_pic2grb(picosCon.ub) 

468 

469 if self.matint: 

470 # Load defining expressions. 

471 gurobiLHS = self._affexp_pic2grb_matint(picosCon.ne) 

472 

473 # Add auxiliary variables for both sides. 

474 gurobiLHSVar = self.int.addMVar(n, lb=-gurobi.GRB.INFINITY) 

475 gurobiRHSVar = self.int.addMVar(1) 

476 

477 # Add constraints to identify auxiliary variables with expressions. 

478 gurobiLHSCon = self.int.addConstr(gurobiLHSVar == gurobiLHS) 

479 gurobiRHSCon = self.int.addConstr(gurobiRHSVar == gurobiRHS) 

480 

481 # Add a quadratic constraint over the auxiliary variables that 

482 # represents the PICOS second order cone constraint itself. 

483 gurobiQuadLHS = gurobiLHSVar.__matmul__(gurobiLHSVar) 

484 gurobiQuadRHS = gurobiRHSVar.__matmul__(gurobiRHSVar) 

485 gurobiQuadCon = self.int.addConstr(gurobiQuadLHS <= gurobiQuadRHS) 

486 

487 # Collect auxiliary objects. 

488 auxCons = [gurobiLHSCon, gurobiRHSCon, gurobiQuadCon] 

489 auxVars = [gurobiLHSVar, gurobiRHSVar] 

490 else: 

491 # Load defining expressions. 

492 gurobiLHS = self._affexp_pic2grb_legacy(picosCon.ne) 

493 

494 # Add auxiliary variables: One for every dimension of the left hand 

495 # side of the PICOS constraint and one for its right hand side. 

496 gurobiLHSVarsDict = self.int.addVars( 

497 n, lb=-gurobi.GRB.INFINITY, ub=gurobi.GRB.INFINITY) 

498 gurobiLHSVars = gurobiLHSVarsDict.values() 

499 gurobiRHSVar = self.int.addVar(lb=0.0, ub=gurobi.GRB.INFINITY) 

500 

501 # Add constraints that identify the left hand side Gurobi auxiliary 

502 # variables with entries of the PICOS left hand side expression. 

503 gurobiLHSDict = dict(enumerate(gurobiLHS)) 

504 gurobiLHSConsDict = self.int.addConstrs( 

505 gurobiLHSVarsDict[d] == gurobiLHSDict[d] for d in range(n)) 

506 gurobiLHSCons = gurobiLHSConsDict.values() 

507 

508 # Add a constraint that identifies the right hand side Gurobi 

509 # auxiliary variable with the PICOS right hand side expression. 

510 gurobiRHSCon = self.int.addLConstr( 

511 gurobiRHSVar, gurobi.GRB.EQUAL, gurobiRHS) 

512 

513 # Add a quadratic constraint over the auxiliary variables that 

514 # represents the PICOS second order cone constraint itself. 

515 quadExpr = gurobi.QuadExpr() 

516 quadExpr.addTerms([1.0] * n, gurobiLHSVars, gurobiLHSVars) 

517 gurobiQuadCon = self.int.addQConstr( 

518 quadExpr, gurobi.GRB.LESS_EQUAL, gurobiRHSVar * gurobiRHSVar) 

519 

520 # Collect auxiliary objects. 

521 auxCons = list(gurobiLHSCons) + [gurobiRHSCon, gurobiQuadCon] 

522 auxVars = list(gurobiLHSVars) + [gurobiRHSVar] 

523 

524 return self.GurobiMetaConstraint(auxCons=auxCons, auxVars=auxVars) 

525 

526 def _import_rscone_constraint(self, picosCon): 

527 import gurobipy as gurobi 

528 

529 assert isinstance(picosCon, RSOCConstraint) 

530 

531 n = len(picosCon.ne) 

532 

533 # Load defining expressions. 

534 gurobiRHS = ( 

535 self._scalar_affexp_pic2grb(picosCon.ub1), 

536 self._scalar_affexp_pic2grb(picosCon.ub2)) 

537 

538 if self.matint: 

539 # Load defining expressions. 

540 gurobiLHS = self._affexp_pic2grb_matint(picosCon.ne) 

541 

542 # Add auxiliary variables for both sides. 

543 gurobiLHSVar = self.int.addMVar(n, lb=-gurobi.GRB.INFINITY) 

544 gurobiRHSVars = (self.int.addMVar(1), self.int.addMVar(1)) 

545 

546 # Add constraints to identify auxiliary variables with expressions. 

547 gurobiLHSCon = self.int.addConstr(gurobiLHSVar == gurobiLHS) 

548 gurobiRHSConsDict = self.int.addConstrs( 

549 gurobiRHSVars[i] == gurobiRHS[i] for i in range(2)) 

550 gurobiRHSCons = gurobiRHSConsDict.values() 

551 

552 # Add a quadratic constraint over the auxiliary variables that 

553 # represents the PICOS rotated second order cone constraint itself. 

554 gurobiQuadLHS = gurobiLHSVar.__matmul__(gurobiLHSVar) 

555 gurobiQuadRHS = gurobiRHSVars[0].__matmul__(gurobiRHSVars[1]) 

556 gurobiQuadCon = self.int.addConstr(gurobiQuadLHS <= gurobiQuadRHS) 

557 

558 # Collect auxiliary objects. 

559 auxCons = [gurobiLHSCon] + list(gurobiRHSCons) + [gurobiQuadCon] 

560 auxVars = [gurobiLHSVar] + list(gurobiRHSVars) 

561 else: 

562 # Load defining expressions. 

563 gurobiLHS = self._affexp_pic2grb_legacy(picosCon.ne) 

564 

565 # Add auxiliary variables: One for every dimension of the left hand 

566 # side of the PICOS constraint and one for its right hand side. 

567 gurobiLHSVarsDict = self.int.addVars( 

568 n, lb=-gurobi.GRB.INFINITY, ub=gurobi.GRB.INFINITY) 

569 gurobiLHSVars = gurobiLHSVarsDict.values() 

570 gurobiRHSVars = self.int.addVars( 

571 2, lb=0.0, ub=gurobi.GRB.INFINITY).values() 

572 

573 # Add constraints that identify the left hand side Gurobi auxiliary 

574 # variables with entries of the PICOS left hand side expression. 

575 gurobiLHSDict = dict(enumerate(gurobiLHS)) 

576 gurobiLHSConsDict = self.int.addConstrs( 

577 gurobiLHSVarsDict[d] == gurobiLHSDict[d] for d in range(n)) 

578 gurobiLHSCons = gurobiLHSConsDict.values() 

579 

580 # Add two constraints that identify the right hand side Gurobi 

581 # auxiliary variables with the PICOS right hand side expressions. 

582 gurobiRHSConsDict = self.int.addConstrs( 

583 gurobiRHSVars[i] == gurobiRHS[i] for i in (0, 1)) 

584 gurobiRHSCons = gurobiRHSConsDict.values() 

585 

586 # Add a quadratic constraint over the auxiliary variables that 

587 # represents the PICOS second order cone constraint itself. 

588 quadExpr = gurobi.QuadExpr() 

589 quadExpr.addTerms([1.0] * n, gurobiLHSVars, gurobiLHSVars) 

590 gurobiQuadCon = self.int.addQConstr(quadExpr, gurobi.GRB.LESS_EQUAL, 

591 gurobiRHSVars[0] * gurobiRHSVars[1]) 

592 

593 # Collect auxiliary objects. 

594 auxCons = ( 

595 list(gurobiLHSCons) + list(gurobiRHSCons) + [gurobiQuadCon]) 

596 auxVars = list(gurobiLHSVars) + list(gurobiRHSVars) 

597 

598 return self.GurobiMetaConstraint(auxCons=auxCons, auxVars=auxVars) 

599 

600 def _import_constraint(self, picosCon): 

601 if isinstance(picosCon, AffineConstraint): 

602 self._gurobiLinCon[picosCon] = \ 

603 self._import_linear_constraint(picosCon) 

604 elif isinstance(picosCon, NonconvexQuadraticConstraint): 

605 self._gurobiQuadCon[picosCon] = \ 

606 self._import_quad_constraint(picosCon) 

607 elif isinstance(picosCon, SOCConstraint): 

608 self._gurobiConicCon[picosCon] = \ 

609 self._import_socone_constraint(picosCon) 

610 elif isinstance(picosCon, RSOCConstraint): 

611 self._gurobiConicCon[picosCon] = \ 

612 self._import_rscone_constraint(picosCon) 

613 else: 

614 assert isinstance(picosCon, DummyConstraint), \ 

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

616 picosCon.__class__.__name__) 

617 

618 def _remove_constraint(self, picosCon): 

619 if isinstance(picosCon, AffineConstraint): 

620 self.int.remove(self._gurobiLinCon.pop(picosCon)) 

621 elif isinstance(picosCon, NonconvexQuadraticConstraint): 

622 self.int.remove(self._gurobiQuadCon.pop(picosCon)) 

623 elif isinstance(picosCon, (SOCConstraint, RSOCConstraint)): 

624 metaCon = self._gurobiConicCon.pop(picosCon) 

625 

626 self.int.remove(metaCon.auxCons) 

627 self.int.remove(metaCon.auxVars) 

628 else: 

629 assert isinstance(picosCon, DummyConstraint), \ 

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

631 picosCon.__class__.__name__) 

632 

633 def _import_objective(self): 

634 import gurobipy as gurobi 

635 

636 picosSense, picosObjective = self.ext.no 

637 

638 # Retrieve objective sense. 

639 if picosSense == "min": 

640 gurobiSense = gurobi.GRB.MINIMIZE 

641 else: 

642 assert picosSense == "max" 

643 gurobiSense = gurobi.GRB.MAXIMIZE 

644 

645 # Retrieve objective function. 

646 if isinstance(picosObjective, AffineExpression): 

647 gurobiObjective = self._scalar_affexp_pic2grb(picosObjective) 

648 else: 

649 assert isinstance(picosObjective, QuadraticExpression) 

650 gurobiObjective = self._quadexp_pic2grb(picosObjective) 

651 

652 self.int.setObjective(gurobiObjective, gurobiSense) 

653 

654 def _import_problem(self): 

655 import gurobipy as gurobi 

656 

657 # Create a problem instance. 

658 if self._license_warnings: 

659 self.int = gurobi.Model() 

660 else: 

661 with self._enforced_verbosity(): 

662 self.int = gurobi.Model() 

663 

664 # Import variables. 

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

666 self._import_variable(variable) 

667 

668 # Import constraints. 

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

670 self._import_constraint(constraint) 

671 

672 # Set objective. 

673 self._import_objective() 

674 

675 def _update_problem(self): 

676 for oldConstraint in self._removed_constraints(): 

677 self._remove_constraint(oldConstraint) 

678 

679 for oldVariable in self._removed_variables(): 

680 self._remove_variable(oldVariable) 

681 

682 for newVariable in self._new_variables(): 

683 self._import_variable(newVariable) 

684 

685 for newConstraint in self._new_constraints(): 

686 self._import_constraint(newConstraint) 

687 

688 if self._objective_has_changed(): 

689 self._import_objective() 

690 

691 def _solve(self): 

692 import gurobipy as gurobi 

693 

694 # Reset options. 

695 # NOTE: OutputFlag = 0 prevents resetParams from printing to console. 

696 self.int.Params.OutputFlag = 0 

697 self.int.resetParams() 

698 

699 # verbosity 

700 self.int.Params.OutputFlag = 1 if self.verbosity() > 0 else 0 

701 

702 # abs_prim_fsb_tol 

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

704 self.int.Params.FeasibilityTol = self.ext.options.abs_prim_fsb_tol 

705 

706 # abs_dual_fsb_tol 

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

708 self.int.Params.OptimalityTol = self.ext.options.abs_dual_fsb_tol 

709 

710 # rel_ipm_opt_tol 

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

712 self.int.Params.BarConvTol = self.ext.options.rel_ipm_opt_tol 

713 

714 # HACK: Work around low precision (conic) quadratic duals. 

715 self.int.Params.BarQCPConvTol = \ 

716 0.01 * self.ext.options.rel_ipm_opt_tol 

717 

718 # abs_bnb_opt_tol 

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

720 self.int.Params.MIPGapAbs = self.ext.options.abs_bnb_opt_tol 

721 

722 # rel_bnb_opt_tol 

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

724 self.int.Params.MIPGap = self.ext.options.rel_bnb_opt_tol 

725 

726 # integrality_tol 

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

728 self.int.Params.IntFeasTol = self.ext.options.integrality_tol 

729 

730 # markowitz_tol 

731 if self.ext.options.markowitz_tol is not None: 

732 self.int.Params.MarkowitzTol = self.ext.options.markowitz_tol 

733 

734 # max_iterations 

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

736 self.int.Params.BarIterLimit = self.ext.options.max_iterations 

737 self.int.Params.IterationLimit = self.ext.options.max_iterations 

738 

739 _lpm = {"interior": 2, "psimplex": 0, "dsimplex": 1} 

740 

741 # lp_node_method 

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

743 value = self.ext.options.lp_node_method 

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

745 self.int.Params.SiftMethod = _lpm[value] 

746 

747 # lp_root_method 

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

749 value = self.ext.options.lp_root_method 

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

751 self.int.Params.Method = _lpm[value] 

752 

753 # timelimit 

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

755 self.int.Params.TimeLimit = self.ext.options.timelimit 

756 

757 # max_fsb_nodes 

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

759 self.int.Params.SolutionLimit = self.ext.options.max_fsb_nodes 

760 

761 # hotstart 

762 if self.ext.options.hotstart: 

763 self._import_variable_values() 

764 else: 

765 self._reset_variable_values() 

766 

767 # Handle Gurobi-specific options. 

768 for key, value in self.ext.options.gurobi_params.items(): 

769 if not self.int.getParamInfo(key): 

770 self._handle_bad_solver_specific_option_key(key) 

771 

772 try: 

773 self.int.setParam(key, value) 

774 except TypeError as error: 

775 self._handle_bad_solver_specific_option_value(key, value, error) 

776 

777 # Handle unsupported options. 

778 self._handle_unsupported_option("treememory") 

779 

780 # Extend functionality for continuous problems. 

781 if self.ext.is_continuous(): 

782 # Compute duals also for QPs and QC(Q)Ps. 

783 if self.ext.options.duals is not False: 

784 self.int.setParam(gurobi.GRB.Param.QCPDual, 1) 

785 

786 # Allow nonconvex quadratic objectives. 

787 # TODO: Allow querying self.ext.objective directly. 

788 # TODO: Check if this should/must be set also for Gurobi >= 9. 

789 if self.ext.footprint.nonconvex_quadratic_objective: 

790 self.int.setParam(gurobi.GRB.Param.NonConvex, 2) 

791 

792 # Attempt to solve the problem. 

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

794 try: 

795 self.int.optimize() 

796 except gurobi.GurobiError as error: 

797 if error.errno == gurobi.GRB.Error.Q_NOT_PSD: 

798 self._handle_continuous_nonconvex_error(error) 

799 else: 

800 raise 

801 

802 # Retrieve primals. 

803 primals = {} 

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

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

806 try: 

807 if self.matint: 

808 value = cvxopt.matrix(self._gurobiVar[picosVar].X) 

809 else: 

810 o = self._gurobiVarOffset[picosVar] 

811 d = picosVar.dim 

812 

813 value = [v.X for v in self._gurobiVars[o:o + d]] 

814 except (AttributeError, gurobi.GurobiError): 

815 # NOTE: AttributeError is raised for gurobipy.Var, 

816 # gurobi.GurobiError for gurobipy.MVar. 

817 primals[picosVar] = None 

818 else: 

819 primals[picosVar] = value 

820 

821 # Retrieve duals. 

822 duals = {} 

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

824 for picosCon in self.ext.constraints.values(): 

825 if isinstance(picosCon, DummyConstraint): 

826 duals[picosCon] = cvxopt.spmatrix([], [], [], picosCon.size) 

827 continue 

828 

829 # HACK: Work around gurobiCon.getAttr(gurobi.GRB.Attr.Pi) 

830 # printing a newline to console when it raises an 

831 # AttributeError and OutputFlag is enabled. This is a 

832 # WONTFIX on Gurobi's end (PICOS #264, Gurobi #14248). 

833 # TODO: Check if this also happens for urobiCon.Pi, which is now 

834 # used for both interfaces. 

835 oldOutput = self.int.Params.OutputFlag 

836 self.int.Params.OutputFlag = 0 

837 

838 try: 

839 if isinstance(picosCon, AffineConstraint): 

840 gurobiCon = self._gurobiLinCon[picosCon] 

841 

842 # HACK: Seee _import_linear_constraint. 

843 if not self.matint or isinstance(gurobiCon, list): 

844 gurobiDual = [c.Pi for c in gurobiCon] 

845 else: 

846 gurobiDual = gurobiCon.Pi 

847 

848 picosDual = cvxopt.matrix(gurobiDual, picosCon.size) 

849 

850 if not picosCon.is_increasing(): 

851 picosDual = -picosDual 

852 elif isinstance(picosCon, SOCConstraint): 

853 gurobiMetaCon = self._gurobiConicCon[picosCon] 

854 

855 if self.matint: 

856 ne, ub, _ = gurobiMetaCon.auxCons 

857 dual = numpy.hstack([ub.Pi, ne.Pi]) 

858 picosDual = cvxopt.matrix(dual) 

859 else: 

860 n = len(picosCon.ne) 

861 assert len(gurobiMetaCon.auxCons) == n + 2 

862 

863 ne = gurobiMetaCon.auxCons[:n] 

864 ub = gurobiMetaCon.auxCons[n] 

865 

866 z, lbd = [c.Pi for c in ne], ub.Pi 

867 picosDual = cvxopt.matrix([lbd] + z) 

868 elif isinstance(picosCon, RSOCConstraint): 

869 gurobiMetaCon = self._gurobiConicCon[picosCon] 

870 

871 if self.matint: 

872 ne, ub1, ub2, _ = gurobiMetaCon.auxCons 

873 dual = numpy.hstack([ub1.Pi, ub2.Pi, ne.Pi]) 

874 picosDual = cvxopt.matrix(dual) 

875 else: 

876 n = len(picosCon.ne) 

877 assert len(gurobiMetaCon.auxCons) == n + 3 

878 

879 ne = gurobiMetaCon.auxCons[:n] 

880 ub1 = gurobiMetaCon.auxCons[n] 

881 ub2 = gurobiMetaCon.auxCons[n + 1] 

882 

883 z, a, b = [c.Pi for c in ne], ub1.Pi, ub2.Pi 

884 picosDual = cvxopt.matrix([a] + [b] + z) 

885 elif isinstance(picosCon, NonconvexQuadraticConstraint): 

886 picosDual = None 

887 else: 

888 assert isinstance(picosCon, DummyConstraint), \ 

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

890 picosCon.__class__.__name__) 

891 

892 # Flip sign based on objective sense. 

893 if picosDual and self.ext.no.direction == "min": 

894 picosDual = -picosDual 

895 except (AttributeError, gurobi.GurobiError): 

896 # NOTE: AttributeError is raised for gurobipy.Constr, 

897 # gurobi.GurobiError for gurobipy.MConstr. 

898 duals[picosCon] = None 

899 else: 

900 duals[picosCon] = picosDual 

901 

902 # HACK: See above. Also: Silence Gurobi while enabling output. 

903 if oldOutput != 0: 

904 with self._enforced_verbosity(noStdOutAt=float("inf")): 

905 self.int.Params.OutputFlag = oldOutput 

906 

907 # Retrieve objective value. 

908 try: 

909 value = self.int.ObjVal 

910 except AttributeError: 

911 value = None 

912 

913 # Retrieve solution status. 

914 statusCode = self.int.Status 

915 if statusCode == gurobi.GRB.Status.LOADED: 

916 raise RuntimeError("Gurobi claims to have just loaded the problem " 

917 "while PICOS expects the solution search to have terminated.") 

918 elif statusCode == gurobi.GRB.Status.OPTIMAL: 

919 primalStatus = SS_OPTIMAL 

920 dualStatus = SS_OPTIMAL 

921 problemStatus = PS_FEASIBLE 

922 elif statusCode == gurobi.GRB.Status.INFEASIBLE: 

923 primalStatus = SS_INFEASIBLE 

924 dualStatus = SS_UNKNOWN 

925 problemStatus = PS_INFEASIBLE 

926 elif statusCode == gurobi.GRB.Status.INF_OR_UNBD: 

927 primalStatus = SS_UNKNOWN 

928 dualStatus = SS_UNKNOWN 

929 problemStatus = PS_INF_OR_UNB 

930 elif statusCode == gurobi.GRB.Status.UNBOUNDED: 

931 primalStatus = SS_UNKNOWN 

932 dualStatus = SS_INFEASIBLE 

933 problemStatus = PS_UNBOUNDED 

934 elif statusCode == gurobi.GRB.Status.CUTOFF: 

935 # "Optimal objective for model was proven to be worse than the value 

936 # specified in the Cutoff parameter. No solution information is 

937 # available." 

938 primalStatus = SS_PREMATURE 

939 dualStatus = SS_PREMATURE 

940 problemStatus = PS_UNKNOWN 

941 elif statusCode == gurobi.GRB.Status.ITERATION_LIMIT: 

942 primalStatus = SS_PREMATURE 

943 dualStatus = SS_PREMATURE 

944 problemStatus = PS_UNKNOWN 

945 elif statusCode == gurobi.GRB.Status.NODE_LIMIT: 

946 primalStatus = SS_PREMATURE 

947 dualStatus = SS_EMPTY # Applies only to mixed integer problems. 

948 problemStatus = PS_UNKNOWN 

949 elif statusCode == gurobi.GRB.Status.TIME_LIMIT: 

950 primalStatus = SS_PREMATURE 

951 dualStatus = SS_PREMATURE 

952 problemStatus = PS_UNKNOWN 

953 elif statusCode == gurobi.GRB.Status.SOLUTION_LIMIT: 

954 primalStatus = SS_PREMATURE 

955 dualStatus = SS_PREMATURE 

956 problemStatus = PS_UNKNOWN 

957 elif statusCode == gurobi.GRB.Status.INTERRUPTED: 

958 primalStatus = SS_PREMATURE 

959 dualStatus = SS_PREMATURE 

960 problemStatus = PS_UNKNOWN 

961 elif statusCode == gurobi.GRB.Status.NUMERIC: 

962 primalStatus = SS_UNKNOWN 

963 dualStatus = SS_UNKNOWN 

964 problemStatus = PS_UNSTABLE 

965 elif statusCode == gurobi.GRB.Status.SUBOPTIMAL: 

966 # "Unable to satisfy optimality tolerances; a sub-optimal solution 

967 # is available." 

968 primalStatus = SS_FEASIBLE 

969 dualStatus = SS_FEASIBLE 

970 problemStatus = PS_FEASIBLE 

971 elif statusCode == gurobi.GRB.Status.INPROGRESS: 

972 raise RuntimeError("Gurobi claims solution search to be 'in " 

973 "progress' while PICOS expects it to have terminated.") 

974 elif statusCode == gurobi.GRB.Status.USER_OBJ_LIMIT: 

975 # "User specified an objective limit (a bound on either the best 

976 # objective or the best bound), and that limit has been reached." 

977 primalStatus = SS_FEASIBLE 

978 dualStatus = SS_EMPTY # Applies only to mixed integer problems. 

979 problemStatus = PS_FEASIBLE 

980 else: 

981 primalStatus = SS_UNKNOWN 

982 dualStatus = SS_UNKNOWN 

983 problemStatus = PS_UNKNOWN 

984 

985 return self._make_solution( 

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

987 

988 

989# -------------------------------------- 

990__all__ = api_end(_API_START, globals())