Coverage for picos/solvers/solver_cplex.py: 79.96%

Shortcuts on this page

r m x   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

514 statements  

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:`CPLEXSolver`.""" 

20 

21import time 

22from collections import namedtuple 

23 

24import cvxopt 

25 

26from ..apidoc import api_end, api_start 

27from ..constraints import (AffineConstraint, ConvexQuadraticConstraint, 

28 DummyConstraint, RSOCConstraint, SOCConstraint) 

29from ..expressions import (CONTINUOUS_VARTYPES, AffineExpression, 

30 BinaryVariable, IntegerVariable, 

31 QuadraticExpression) 

32from ..modeling.footprint import Specification 

33from ..modeling.solution import (PS_FEASIBLE, PS_ILLPOSED, PS_INF_OR_UNB, 

34 PS_INFEASIBLE, PS_UNBOUNDED, PS_UNKNOWN, 

35 PS_UNSTABLE, SS_EMPTY, SS_FAILURE, 

36 SS_FEASIBLE, SS_INFEASIBLE, SS_OPTIMAL, 

37 SS_PREMATURE, SS_UNKNOWN) 

38from .solver import (ConflictingOptionsError, DependentOptionError, Solver, 

39 UnsupportedOptionError) 

40 

41_API_START = api_start(globals()) 

42# ------------------------------- 

43 

44 

45#: Maps CPLEX status code to PICOS status triples. 

46CPLEX_STATUS_CODES = { 

47 # primal status, dual status, problem status 

481: (SS_OPTIMAL, SS_OPTIMAL, PS_FEASIBLE), # CPX_STAT_OPTIMAL 

492: (SS_UNKNOWN, SS_INFEASIBLE, PS_UNBOUNDED), # CPX_STAT_UNBOUNDED 

503: (SS_INFEASIBLE, SS_UNKNOWN, PS_INFEASIBLE), # CPX_STAT_INFEASIBLE 

514: (SS_UNKNOWN, SS_UNKNOWN, PS_INF_OR_UNB), # CPX_STAT_INForUNBD 

525: (SS_INFEASIBLE, SS_UNKNOWN, PS_UNSTABLE), # CPX_STAT_OPTIMAL_INFEAS 

536: (SS_UNKNOWN, SS_UNKNOWN, PS_UNSTABLE), # CPX_STAT_NUM_BEST 

54 # 7—9 are not defined. 

5510: (SS_PREMATURE, SS_PREMATURE, PS_UNKNOWN), # CPX_STAT_ABORT_IT_LIM 

5611: (SS_PREMATURE, SS_PREMATURE, PS_UNKNOWN), # CPX_STAT_ABORT_TIME_LIM 

5712: (SS_PREMATURE, SS_PREMATURE, PS_UNKNOWN), # CPX_STAT_ABORT_OBJ_LIM 

5813: (SS_PREMATURE, SS_PREMATURE, PS_UNKNOWN), # CPX_STAT_ABORT_USER 

59 # 14—19 seem irrelevant (CPX_STAT_*_RELAXED_*). 

6020: (SS_UNKNOWN, SS_UNKNOWN, PS_ILLPOSED), # …_OPTIMAL_FACE_UNBOUNDED 

6121: (SS_PREMATURE, SS_PREMATURE, PS_UNKNOWN), # …_ABORT_PRIM_OBJ_LIM 

6222: (SS_PREMATURE, SS_PREMATURE, PS_UNKNOWN), # …_ABORT_DUAL_OBJ_LIM 

6323: (SS_FEASIBLE, SS_FEASIBLE, PS_FEASIBLE), # CPX_STAT_FEASIBLE 

64 # 24 irrelevant (CPX_STAT_FIRSTORDER). 

6525: (SS_PREMATURE, SS_PREMATURE, PS_UNKNOWN), # …_ABORT_DETTIME_LIM 

66 # 26—29 are not defined. 

67 # 30—39 seem irrelevant (CPX_STAT_CONFLICT_*). 

68 # 40—100 are not defined. 

69101: (SS_OPTIMAL, SS_EMPTY, PS_FEASIBLE), # CPXMIP_OPTIMAL 

70102: (SS_OPTIMAL, SS_EMPTY, PS_FEASIBLE), # CPXMIP_OPTIMAL_TOL 

71103: (SS_INFEASIBLE, SS_EMPTY, PS_INFEASIBLE), # CPXMIP_INFEASIBLE 

72104: (SS_PREMATURE, SS_EMPTY, PS_UNKNOWN), # CPXMIP_SOL_LIM ? 

73105: (SS_FEASIBLE, SS_EMPTY, PS_FEASIBLE), # CPXMIP_NODE_LIM_FEAS 

74106: (SS_PREMATURE, SS_EMPTY, PS_UNKNOWN), # CPXMIP_NODE_LIM_INFEAS 

75107: (SS_FEASIBLE, SS_EMPTY, PS_FEASIBLE), # CPXMIP_TIME_LIM_FEAS 

76108: (SS_PREMATURE, SS_EMPTY, PS_UNKNOWN), # CPXMIP_TIME_LIM_INFEAS 

77109: (SS_FEASIBLE, SS_EMPTY, PS_FEASIBLE), # CPXMIP_FAIL_FEAS 

78110: (SS_FAILURE, SS_EMPTY, PS_UNKNOWN), # CPXMIP_FAIL_INFEAS 

79111: (SS_FEASIBLE, SS_EMPTY, PS_FEASIBLE), # CPXMIP_MEM_LIM_FEAS 

80112: (SS_PREMATURE, SS_EMPTY, PS_UNKNOWN), # CPXMIP_MEM_LIM_INFEAS 

81113: (SS_FEASIBLE, SS_EMPTY, PS_FEASIBLE), # CPXMIP_ABORT_FEAS 

82114: (SS_PREMATURE, SS_EMPTY, PS_UNKNOWN), # CPXMIP_ABORT_INFEAS 

83115: (SS_INFEASIBLE, SS_EMPTY, PS_UNSTABLE), # CPXMIP_OPTIMAL_INFEAS 

84116: (SS_FEASIBLE, SS_EMPTY, PS_FEASIBLE), # CPXMIP_FAIL_FEAS_NO_TREE 

85117: (SS_FAILURE, SS_EMPTY, PS_UNKNOWN), # …_FAIL_INFEAS_NO_TREE 

86118: (SS_UNKNOWN, SS_EMPTY, PS_UNBOUNDED), # CPXMIP_UNBOUNDED 

87119: (SS_UNKNOWN, SS_EMPTY, PS_INF_OR_UNB), # CPXMIP_INForUNBD 

88 # 120—126 seem irrelevant (CPXMIP_*_RELAXED_*). 

89127: (SS_FEASIBLE, SS_EMPTY, PS_FEASIBLE), # CPXMIP_FEASIBLE 

90128: (SS_OPTIMAL, SS_EMPTY, PS_FEASIBLE), # …_POPULATESOL_LIM ? 

91129: (SS_OPTIMAL, SS_EMPTY, PS_FEASIBLE), # …_OPTIMAL_POPULATED ? 

92130: (SS_OPTIMAL, SS_EMPTY, PS_FEASIBLE), # …_OPTIMAL_POPULATED_TOL ? 

93131: (SS_FEASIBLE, SS_EMPTY, PS_FEASIBLE), # CPXMIP_DETTIME_LIM_FEAS 

94132: (SS_PREMATURE, SS_EMPTY, PS_UNKNOWN), # CPXMIP_DETTIME_LIM_INFEAS 

95} 

96 

97 

98class CPLEXSolver(Solver): 

99 """Interface to the CPLEX solver via its official Python interface. 

100 

101 .. note :: 

102 Names are used instead of indices for identifying both variables and 

103 constraints since indices can change if the CPLEX instance is modified. 

104 """ 

105 

106 # NOTE: When making changes, also see the section in _solve that tells CPLEX 

107 # the problem type. 

108 SUPPORTED = Specification( 

109 objectives=[ 

110 AffineExpression, 

111 QuadraticExpression], 

112 constraints=[ 

113 DummyConstraint, 

114 AffineConstraint, 

115 SOCConstraint, 

116 RSOCConstraint, 

117 ConvexQuadraticConstraint]) 

118 

119 NONCONVEX_QP = Specification( 

120 objectives=[QuadraticExpression], 

121 constraints=[DummyConstraint, AffineConstraint]) 

122 

123 MetaConstraint = namedtuple("MetaConstraint", ("con", "dim")) 

124 

125 @classmethod 

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

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

128 result = Solver.supports(footprint, explain) 

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

130 return result 

131 

132 # Support QPs and MIQPs with a nonconvex objective. 

133 # NOTE: SUPPORTED fully excludes nonconvex quadratic constraints. This 

134 # further excludes QCQPs and MIQCQPs with a nonconvex objective. 

135 # TODO: See which of the excluded cases can be supported as well. 

136 if footprint.nonconvex_quadratic_objective \ 

137 and footprint not in cls.NONCONVEX_QP: 

138 if explain: 

139 return (False, "(MI)QCQPs with nonconvex objective.") 

140 else: 

141 return False 

142 

143 if footprint not in cls.SUPPORTED: 

144 if explain: 

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

146 else: 

147 return False 

148 

149 return (True, None) if explain else True 

150 

151 @classmethod 

152 def default_penalty(cls): 

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

154 return 0.0 # Commercial solver. 

155 

156 @classmethod 

157 def test_availability(cls): 

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

159 cls.check_import("cplex") 

160 

161 @classmethod 

162 def names(cls): 

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

164 return "cplex", "CPLEX", "IBM ILOG CPLEX Optimization Studio", None 

165 

166 @classmethod 

167 def is_free(cls): 

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

169 return False 

170 

171 def __init__(self, problem): 

172 """Initialize a CPLEX solver interface. 

173 

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

175 """ 

176 super(CPLEXSolver, self).__init__(problem) 

177 

178 self._cplexVar = dict(start=dict(), length=dict()) 

179 """Maps a PICOS variable to a CPLEX start index and length.""" 

180 

181 self._cplexLinCon = dict(start=dict(), length=dict()) 

182 """Maps a PICOS linear constraint to a CPLEX start index and length.""" 

183 

184 self._cplexQuadCon = dict(start=dict(), length=dict()) 

185 """Maps a PICOS quadr. constraint to a CPLEX start index and length.""" 

186 

187 self._cplexMetaCon = dict() 

188 """Maps PICOS (rotated) second order conic constraints to a named tuple. 

189 

190 The tuple has a ``dim`` property containing the linear auxiliary 

191 variable dimension and can thus be used as a key of ``self._cplexVar``. 

192 """ 

193 

194 def __del__(self): 

195 if self.int is not None: 

196 self.int.end() 

197 

198 def reset_problem(self): 

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

200 if self.int is not None: 

201 self.int.end() 

202 self.int = None 

203 

204 self._cplexVar["start"].clear() 

205 self._cplexVar["length"].clear() 

206 

207 self._cplexLinCon["start"].clear() 

208 self._cplexLinCon["length"].clear() 

209 

210 self._cplexQuadCon["start"].clear() 

211 self._cplexQuadCon["length"].clear() 

212 

213 self._cplexMetaCon.clear() 

214 

215 @classmethod 

216 def _register(cls, registry, key, indices): 

217 if isinstance(indices, int): 

218 start, length = indices, 1 

219 else: 

220 start, length = indices[0], len(indices) 

221 

222 # Expect that indices are consecutive. 

223 assert isinstance(indices, range) \ 

224 or tuple(indices) == tuple(range(start, start + length)), \ 

225 "Not consecutive: {}".format(indices) 

226 

227 registry["start"][key] = start 

228 registry["length"][key] = length 

229 

230 @classmethod 

231 def _lookup(cls, registry, key): 

232 start = registry["start"][key] 

233 length = registry["length"][key] 

234 

235 return list(range(start, start + length)) 

236 

237 @classmethod 

238 def _unregister(cls, registry, key): 

239 starts = registry["start"] 

240 

241 start = starts.pop(key) 

242 length = registry["length"].pop(key) 

243 indices = list(range(start, start + length)) 

244 

245 for other in starts: 

246 if starts[other] > start: 

247 starts[other] -= length 

248 

249 return indices 

250 

251 def _import_variable(self, picosVar): 

252 import cplex 

253 

254 dim = picosVar.dim 

255 

256 # Retrieve types. 

257 if isinstance(picosVar, CONTINUOUS_VARTYPES): 

258 types = dim * self.int.variables.type.continuous 

259 elif isinstance(picosVar, IntegerVariable): 

260 types = dim * self.int.variables.type.integer 

261 elif isinstance(picosVar, BinaryVariable): 

262 types = dim * self.int.variables.type.binary 

263 else: 

264 assert False, "Unexpected variable type." 

265 

266 # Retrieve bounds. 

267 lowerBounds = [-cplex.infinity]*dim 

268 upperBounds = [cplex.infinity]*dim 

269 lower, upper = picosVar.bound_dicts 

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

271 lowerBounds[i] = b 

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

273 upperBounds[i] = b 

274 

275 # Import the variable. 

276 cplexIndices = self.int.variables.add( 

277 lb=lowerBounds, ub=upperBounds, types=types) 

278 

279 # Register the variable. 

280 self._register(self._cplexVar, picosVar, cplexIndices) 

281 

282 def _remove_variable(self, picosVar): 

283 # Unregister the variable. 

284 cplexIndices = self._unregister(self._cplexVar, picosVar) 

285 

286 # Remove the variable. 

287 self.int.variables.delete(cplexIndices) 

288 

289 def _affinexp_pic2cpl(self, picosExpression): 

290 import cplex 

291 

292 for I, V, c in picosExpression.sparse_rows(self._cplexVar["start"]): 

293 yield cplex.SparsePair(ind=I, val=V), c 

294 

295 def _scalar_affinexp_pic2cpl(self, picosExpression): 

296 assert len(picosExpression) == 1 

297 

298 return next(self._affinexp_pic2cpl(picosExpression)) 

299 

300 def _quadexp_pic2cpl(self, picosExpression): 

301 import cplex 

302 

303 assert isinstance(picosExpression, QuadraticExpression) 

304 

305 start = self._cplexVar["start"] 

306 

307 cplexI, cplexJ, cplexV = [], [], [] 

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

309 cplexI.extend(Q.I + start[x]) 

310 cplexJ.extend(Q.J + start[y]) 

311 cplexV.extend(Q.V) 

312 

313 return cplex.SparseTriple(ind1=cplexI, ind2=cplexJ, val=cplexV) 

314 

315 def _import_linear_constraint(self, picosCon): 

316 assert isinstance(picosCon, AffineConstraint) 

317 

318 length = len(picosCon) 

319 

320 # Retrieve left hand side and right hand side expressions. 

321 cplexLHS, cplexRHS = [], [] 

322 for linear, constant in self._affinexp_pic2cpl(picosCon.lmr): 

323 cplexLHS.append(linear) 

324 cplexRHS.append(-constant) 

325 

326 # Retrieve senses. 

327 if picosCon.is_increasing(): 

328 senses = length * "L" 

329 elif picosCon.is_decreasing(): 

330 senses = length * "G" 

331 elif picosCon.is_equality(): 

332 senses = length * "E" 

333 else: 

334 assert False, "Unexpected constraint relation." 

335 

336 # Import the constraint. 

337 cplexIndices = self.int.linear_constraints.add( 

338 lin_expr=cplexLHS, senses=senses, rhs=cplexRHS) 

339 

340 # Register the constraint. 

341 self._register(self._cplexLinCon, picosCon, cplexIndices) 

342 

343 def _import_quad_constraint(self, picosCon): 

344 assert isinstance(picosCon, ConvexQuadraticConstraint) 

345 

346 # Retrieve the affine term. 

347 cplexLinear, cplexRHS = self._scalar_affinexp_pic2cpl(picosCon.le0.aff) 

348 cplexRHS = -cplexRHS 

349 

350 # Retrieve the quadratic term. 

351 cplexQuad = self._quadexp_pic2cpl(picosCon.le0) 

352 

353 # Import the constraint. 

354 cplexIndices = self.int.quadratic_constraints.add( 

355 lin_expr=cplexLinear, quad_expr=cplexQuad, sense="L", rhs=cplexRHS) 

356 

357 # Register the constraint. 

358 self._register(self._cplexQuadCon, picosCon, cplexIndices) 

359 

360 # TODO: Handle SOC → Quadratic via a reformulation. 

361 def _import_socone_constraint(self, picosCon): 

362 import cplex 

363 

364 assert isinstance(picosCon, SOCConstraint) 

365 

366 picosLHS = picosCon.ne 

367 picosRHS = picosCon.ub 

368 picosLHSLen = len(picosLHS) 

369 

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

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

372 cplexRHSVar = self.int.variables.add( 

373 lb=[0.0], ub=[+cplex.infinity], 

374 types=self.int.variables.type.continuous)[0] 

375 cplexLHSVars = self.int.variables.add( 

376 lb=[-cplex.infinity] * picosLHSLen, 

377 ub=[+cplex.infinity] * picosLHSLen, 

378 types=self.int.variables.type.continuous * picosLHSLen) 

379 

380 # Add a constraint that identifies the right hand side CPLEX auxiliary 

381 # variable with the PICOS right hand side scalar expression. 

382 # NOTE: Order (RHS first) matters for dual retrieval. 

383 cplexRHSConLHS, cplexRHSConRHS = \ 

384 self._scalar_affinexp_pic2cpl(-picosRHS) 

385 cplexRHSConRHS = -cplexRHSConRHS 

386 cplexRHSConLHS.ind.append(cplexRHSVar) 

387 cplexRHSConLHS.val.append(1.0) 

388 cplexRHSCon = self.int.linear_constraints.add( 

389 lin_expr=[cplexRHSConLHS], senses="E", rhs=[cplexRHSConRHS])[0] 

390 

391 # Add constraints that identify the left hand side CPLEX auxiliary 

392 # variables with their slice of the PICOS left hand side expression. 

393 # TODO: Possible to get rid of the loop? 

394 cplexLHSConsLHSs, cplexLHSConsRHSs = [], [] 

395 for localConIndex, (localLinExp, localConstant) in \ 

396 enumerate(self._affinexp_pic2cpl(picosLHS)): 

397 localConstant = -localConstant 

398 localLinExp.ind.append(cplexLHSVars[localConIndex]) 

399 localLinExp.val.append(-1.0) 

400 cplexLHSConsLHSs.append(localLinExp) 

401 cplexLHSConsRHSs.append(localConstant) 

402 cplexLHSCons = self.int.linear_constraints.add( 

403 lin_expr=cplexLHSConsLHSs, senses="E" * picosLHSLen, 

404 rhs=cplexLHSConsRHSs) 

405 

406 # Add a quadratic constraint over the auxiliary variables that 

407 # represents the PICOS second order cone constraint itself. 

408 quadIndices = [cplexRHSVar] + list(cplexLHSVars) 

409 quadExpr = cplex.SparseTriple( 

410 ind1=quadIndices, ind2=quadIndices, val=[-1.0] + [1.0]*picosLHSLen) 

411 cplexQuadCon = self.int.quadratic_constraints.add( 

412 quad_expr=quadExpr, sense="L", rhs=0.0) 

413 

414 # Register all auxiliary variables and constraints. 

415 cplexVars = [cplexRHSVar] + list(cplexLHSVars) 

416 cplexLinCons = [cplexRHSCon] + list(cplexLHSCons) 

417 metaCon = self.MetaConstraint(con=picosCon, dim=len(cplexVars)) 

418 

419 self._cplexMetaCon[picosCon] = metaCon 

420 self._register(self._cplexVar, metaCon, cplexVars) 

421 self._register(self._cplexLinCon, metaCon, cplexLinCons) 

422 self._register(self._cplexQuadCon, metaCon, cplexQuadCon) 

423 

424 # TODO: Handle RSOC → Quadratic via a reformulation. 

425 def _import_rscone_constraint(self, picosCon): 

426 import cplex 

427 

428 assert isinstance(picosCon, RSOCConstraint) 

429 

430 picosLHS = picosCon.ne 

431 picosRHS1 = picosCon.ub1 

432 picosRHS2 = picosCon.ub2 

433 picosLHSLen = len(picosLHS) 

434 

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

436 # of the PICOS constraint and two for its right hand side. 

437 cplexRHSVars = self.int.variables.add( 

438 lb=[0.0, 0.0], ub=[+cplex.infinity] * 2, 

439 types=self.int.variables.type.continuous * 2) 

440 cplexLHSVars = self.int.variables.add( 

441 lb=[-cplex.infinity] * picosLHSLen, 

442 ub=[+cplex.infinity] * picosLHSLen, 

443 types=self.int.variables.type.continuous * picosLHSLen) 

444 

445 # Add two constraints that identify the right hand side CPLEX auxiliary 

446 # variables with the PICOS right hand side scalar expressions. 

447 # NOTE: Order (RHS first) matters for dual retrieval. 

448 cplexRHSConsLHSs, cplexRHSConsRHSs = [], [] 

449 for picosRHS, cplexRHSVar in zip((picosRHS1, picosRHS2), cplexRHSVars): 

450 linExp, constant = self._scalar_affinexp_pic2cpl(-picosRHS) 

451 linExp.ind.append(cplexRHSVar) 

452 linExp.val.append(1.0) 

453 constant = -constant 

454 cplexRHSConsLHSs.append(linExp) 

455 cplexRHSConsRHSs.append(constant) 

456 cplexRHSCons = self.int.linear_constraints.add( 

457 lin_expr=cplexRHSConsLHSs, senses="E" * 2, rhs=cplexRHSConsRHSs) 

458 

459 # Add constraints that identify the left hand side CPLEX auxiliary 

460 # variables with their slice of the PICOS left hand side expression. 

461 # TODO: Possible to get rid of the loop? 

462 cplexLHSConsLHSs, cplexLHSConsRHSs = [], [] 

463 for localConIndex, (localLinExp, localConstant) in \ 

464 enumerate(self._affinexp_pic2cpl(picosLHS)): 

465 localLinExp.ind.append(cplexLHSVars[localConIndex]) 

466 localLinExp.val.append(-1.0) 

467 localConstant = -localConstant 

468 cplexLHSConsLHSs.append(localLinExp) 

469 cplexLHSConsRHSs.append(localConstant) 

470 cplexLHSCons = self.int.linear_constraints.add( 

471 lin_expr=cplexLHSConsLHSs, senses="E" * picosLHSLen, 

472 rhs=cplexLHSConsRHSs) 

473 

474 # Add a quadratic constraint over the auxiliary variables that 

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

476 quadExpr = cplex.SparseTriple( 

477 ind1=[cplexRHSVars[0]] + list(cplexLHSVars), 

478 ind2=[cplexRHSVars[1]] + list(cplexLHSVars), 

479 val=[-1.0] + [1.0] * picosLHSLen) 

480 cplexQuadCon = self.int.quadratic_constraints.add( 

481 quad_expr=quadExpr, sense="L", rhs=0.0) 

482 

483 # Register all auxiliary variables and constraints. 

484 cplexVars = list(cplexRHSVars) + list(cplexLHSVars) 

485 cplexLinCons = list(cplexRHSCons) + list(cplexLHSCons) 

486 

487 metaCon = self.MetaConstraint(con=picosCon, dim=len(cplexVars)) 

488 

489 self._cplexMetaCon[picosCon] = metaCon 

490 self._register(self._cplexVar, metaCon, cplexVars) 

491 self._register(self._cplexLinCon, metaCon, cplexLinCons) 

492 self._register(self._cplexQuadCon, metaCon, cplexQuadCon) 

493 

494 def _import_constraint(self, picosCon): 

495 if isinstance(picosCon, AffineConstraint): 

496 self._import_linear_constraint(picosCon) 

497 elif isinstance(picosCon, ConvexQuadraticConstraint): 

498 self._import_quad_constraint(picosCon) 

499 elif isinstance(picosCon, SOCConstraint): 

500 self._import_socone_constraint(picosCon) 

501 elif isinstance(picosCon, RSOCConstraint): 

502 self._import_rscone_constraint(picosCon) 

503 else: 

504 assert isinstance(picosCon, DummyConstraint), \ 

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

506 picosCon.__class__.__name__) 

507 

508 def _remove_constraint(self, picosCon): 

509 if isinstance(picosCon, AffineConstraint): 

510 cplexIndices = self._unregister(self._cplexLinCon, picosCon) 

511 

512 self.int.linear_constraints.delete(cplexIndices) 

513 elif isinstance(picosCon, ConvexQuadraticConstraint): 

514 cplexIndices = self._unregister(self._cplexQuadCon, picosCon) 

515 

516 self.int.quadratic_constraints.delete(cplexIndices) 

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

518 metaCon = self._cplexMetaCon.pop(picosCon) 

519 

520 cplexLinConIndices = self._unregister(self._cplexLinCon, metaCon) 

521 cplexQuadConIndices = self._unregister(self._cplexQuadCon, metaCon) 

522 cplexVarIndices = self._unregister(self._cplexVar, metaCon) 

523 

524 self.int.linear_constraints.delete(cplexLinConIndices) 

525 self.int.quadratic_constraints.delete(cplexQuadConIndices) 

526 self.int.variables.delete(cplexVarIndices) 

527 else: 

528 assert isinstance(picosCon, DummyConstraint), \ 

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

530 picosCon.__class__.__name__) 

531 

532 def _import_affine_objective(self, picosExpression): 

533 assert isinstance(picosExpression, AffineExpression) 

534 assert picosExpression.scalar 

535 

536 # Import constant part. 

537 self.int.objective.set_offset(picosExpression._constant_coef[0]) 

538 

539 # Import linear part. 

540 cplexLinear = [] 

541 for picosVar, coefs in picosExpression._sparse_linear_coefs.items(): 

542 cplexIndices = coefs.J + self._cplexVar["start"][picosVar] 

543 cplexCoefs = list(coefs) 

544 

545 cplexLinear.extend(zip(cplexIndices, cplexCoefs)) 

546 

547 if cplexLinear: 

548 self.int.objective.set_linear(cplexLinear) 

549 

550 def _reset_affine_objective(self): 

551 # Clear constant part. 

552 self.int.objective.set_offset(0.0) 

553 

554 # Clear linear part. 

555 linear = self.int.objective.get_linear() 

556 if any(linear): 

557 self.int.objective.set_linear([(cplexVarIndex, 0.0) 

558 for cplexVarIndex, coef in enumerate(linear) if coef]) 

559 

560 def _import_quadratic_objective(self, picosExpression): 

561 assert isinstance(picosExpression, QuadraticExpression) 

562 

563 # Import affine part of objective function. 

564 self._import_affine_objective(picosExpression.aff) 

565 

566 # Import quadratic part of objective function. 

567 cplexQuadExpression = self._quadexp_pic2cpl(picosExpression) 

568 cplexQuadCoefs = zip( 

569 cplexQuadExpression.ind1, cplexQuadExpression.ind2, 

570 [2.0 * coef for coef in cplexQuadExpression.val]) 

571 self.int.objective.set_quadratic_coefficients(cplexQuadCoefs) 

572 

573 def _reset_quadratic_objective(self): 

574 quadratics = self.int.objective.get_quadratic() 

575 if quadratics: 

576 self.int.objective.set_quadratic( 

577 [(sparsePair.ind, [0]*len(sparsePair.ind)) 

578 for sparsePair in quadratics]) 

579 

580 def _import_objective(self): 

581 picosSense, picosObjective = self.ext.no 

582 

583 # Import objective sense. 

584 if picosSense == "min": 

585 cplexSense = self.int.objective.sense.minimize 

586 else: 

587 assert picosSense == "max" 

588 cplexSense = self.int.objective.sense.maximize 

589 self.int.objective.set_sense(cplexSense) 

590 

591 # Import objective function. 

592 if isinstance(picosObjective, AffineExpression): 

593 self._import_affine_objective(picosObjective) 

594 else: 

595 assert isinstance(picosObjective, QuadraticExpression) 

596 self._import_quadratic_objective(picosObjective) 

597 

598 def _reset_objective(self): 

599 self._reset_affine_objective() 

600 self._reset_quadratic_objective() 

601 

602 def _import_problem(self): 

603 import cplex 

604 

605 # Create a problem instance. 

606 self.int = cplex.Cplex() 

607 

608 # Import variables. 

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

610 self._import_variable(variable) 

611 

612 # Import constraints. 

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

614 self._import_constraint(constraint) 

615 

616 # Set objective. 

617 self._import_objective() 

618 

619 def _update_problem(self): 

620 for oldConstraint in self._removed_constraints(): 

621 self._remove_constraint(oldConstraint) 

622 

623 for oldVariable in self._removed_variables(): 

624 self._remove_variable(oldVariable) 

625 

626 for newVariable in self._new_variables(): 

627 self._import_variable(newVariable) 

628 

629 for newConstraint in self._new_constraints(): 

630 self._import_constraint(newConstraint) 

631 

632 if self._objective_has_changed(): 

633 self._reset_objective() 

634 self._import_objective() 

635 

636 def _solve(self): 

637 import cplex 

638 

639 # Reset options. 

640 self.int.parameters.reset() 

641 

642 o = self.ext.options 

643 p = self.int.parameters 

644 

645 continuous = self.ext.is_continuous() 

646 

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

648 nonconvex_quad_obj = self.ext.footprint.nonconvex_quadratic_objective 

649 

650 # verbosity 

651 verbosity = self.verbosity() 

652 if verbosity <= 0: 

653 # Note that this behaviour disables warning even with a verbosity of 

654 # zero but this is still better than having verbose output for every 

655 # option that is set. 

656 self.int.set_results_stream(None) 

657 else: 

658 p.barrier.display.set(min(2, verbosity)) 

659 p.conflict.display.set(min(2, verbosity)) 

660 p.mip.display.set(min(5, verbosity)) 

661 p.sifting.display.set(min(2, verbosity)) 

662 p.simplex.display.set(min(2, verbosity)) 

663 p.tune.display.set(min(3, verbosity)) 

664 

665 self.int.set_error_stream(None) # Already handled as exceptions. 

666 

667 # abs_prim_fsb_tol 

668 if o.abs_prim_fsb_tol is not None: 

669 p.simplex.tolerances.feasibility.set(o.abs_prim_fsb_tol) 

670 

671 # abs_dual_fsb_tol 

672 if o.abs_dual_fsb_tol is not None: 

673 p.simplex.tolerances.optimality.set(o.abs_dual_fsb_tol) 

674 

675 # rel_prim_fsb_tol, rel_dual_fsb_tol, rel_ipm_opt_tol 

676 convergenceTols = [tol for tol in (o.rel_prim_fsb_tol, 

677 o.rel_dual_fsb_tol, o.rel_ipm_opt_tol) if tol is not None] 

678 if convergenceTols: 

679 convergenceTol = min(convergenceTols) 

680 p.barrier.convergetol.set(convergenceTol) 

681 p.barrier.qcpconvergetol.set(convergenceTol) 

682 

683 # abs_bnb_opt_tol 

684 if o.abs_bnb_opt_tol is not None: 

685 p.mip.tolerances.absmipgap.set(o.abs_bnb_opt_tol) 

686 

687 # rel_bnb_opt_tol 

688 if o.rel_bnb_opt_tol is not None: 

689 p.mip.tolerances.mipgap.set(o.rel_bnb_opt_tol) 

690 

691 # integrality_tol 

692 if o.integrality_tol is not None: 

693 p.mip.tolerances.integrality.set(o.integrality_tol) 

694 

695 # markowitz_tol 

696 if o.markowitz_tol is not None: 

697 p.simplex.tolerances.markowitz.set(o.markowitz_tol) 

698 

699 # max_iterations 

700 if o.max_iterations is not None: 

701 maxit = o.max_iterations 

702 p.barrier.limits.iteration.set(maxit) 

703 p.simplex.limits.iterations.set(maxit) 

704 

705 _lpm = {"interior": 4, "psimplex": 1, "dsimplex": 2} 

706 

707 # lp_node_method 

708 if o.lp_node_method is not None: 

709 assert o.lp_node_method in _lpm, "Unexpected lp_node_method value." 

710 p.mip.strategy.subalgorithm.set(_lpm[o.lp_node_method]) 

711 

712 # lp_root_method 

713 if o.lp_root_method is not None: 

714 assert o.lp_root_method in _lpm, "Unexpected lp_root_method value." 

715 p.lpmethod.set(_lpm[o.lp_root_method]) 

716 

717 # timelimit 

718 if o.timelimit is not None: 

719 p.timelimit.set(o.timelimit) 

720 

721 # treememory 

722 if o.treememory is not None: 

723 p.mip.limits.treememory.set(o.treememory) 

724 

725 # Handle option conflict between "max_fsb_nodes" and "pool_size". 

726 if o.max_fsb_nodes is not None \ 

727 and o.pool_size is not None: 

728 raise ConflictingOptionsError("The options 'max_fsb_nodes' and " 

729 "'pool_size' cannot be used in conjunction.") 

730 

731 # max_fsb_nodes 

732 if o.max_fsb_nodes is not None: 

733 p.mip.limits.solutions.set(o.max_fsb_nodes) 

734 

735 # pool_size 

736 if o.pool_size is not None: 

737 if continuous: 

738 raise UnsupportedOptionError("The option 'pool_size' can only " 

739 "be used with mixed integer problems.") 

740 maxNumSolutions = max(1, int(o.pool_size)) 

741 p.mip.limits.populate.set(maxNumSolutions) 

742 else: 

743 maxNumSolutions = 1 

744 

745 # pool_relgap 

746 if o.pool_rel_gap is not None: 

747 if o.pool_size is None: 

748 raise DependentOptionError("The option 'pool_rel_gap' requires " 

749 "the option 'pool_size'.") 

750 p.mip.pool.relgap.set(o.pool_rel_gap) 

751 

752 # pool_abs_gap 

753 if o.pool_abs_gap is not None: 

754 if o.pool_size is None: 

755 raise DependentOptionError("The option 'pool_abs_gap' requires " 

756 "the option 'pool_size'.") 

757 p.mip.pool.absgap.set(o.pool_abs_gap) 

758 

759 # hotstart 

760 if o.hotstart: 

761 indices, values = [], [] 

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

763 if picosVar.valued: 

764 indices.extend(self._lookup(self._cplexVar, picosVar)) 

765 values.extend(cvxopt.matrix(picosVar.internal_value)) 

766 

767 if indices: 

768 self.int.MIP_starts.add( 

769 cplex.SparsePair(ind=indices, val=values), 

770 self.int.MIP_starts.effort_level.repair) 

771 

772 # Set the optimality target now so that cplex_params may overwrite it. 

773 # This allows solving QPs and MIQPs with a nonconvex objective. 

774 if nonconvex_quad_obj: 

775 p.optimalitytarget.set(3) 

776 

777 # Load a virtual machine config. 

778 if self.ext.options.cplex_vmconfig: 

779 self.int.copy_vmconfig(self.ext.options.cplex_vmconfig) 

780 

781 # Handle CPLEX-specific options. 

782 for key, value in o.cplex_params.items(): 

783 try: 

784 parameter = getattr(self.int.parameters, key) 

785 except AttributeError as error: 

786 self._handle_bad_solver_specific_option_key(key, error) 

787 

788 try: 

789 parameter.set(value) 

790 except cplex.exceptions.errors.CplexError as error: 

791 self._handle_bad_solver_specific_option_value(key, value, error) 

792 

793 # Handle options "cplex_upr_bnd_limit", "cplex_lwr_bnd_limit" and 

794 # "cplex_bnd_monitor" via a CPLEX callback handler. 

795 callback = None 

796 if o.cplex_upr_bnd_limit or o.cplex_lwr_bnd_limit \ 

797 or o.cplex_bnd_monitor: 

798 from cplex.callbacks import MIPInfoCallback 

799 

800 class PicosInfoCallback(MIPInfoCallback): 

801 def __call__(self): 

802 v1 = self.get_incumbent_objective_value() 

803 v2 = self.get_best_objective_value() 

804 ub = max(v1, v2) 

805 lb = min(v1, v2) 

806 if self.bounds is not None: 

807 elapsedTime = time.time() - self.startTime 

808 self.bounds.append((elapsedTime, lb, ub)) 

809 if self.lbound is not None and lb >= self.lbound: 

810 self.printer("The specified lower bound was reached, " 

811 "so PICOS will ask CPLEX to stop the search.") 

812 self.abort() 

813 if self.ubound is not None and ub <= self.ubound: 

814 self.printer("The specified upper bound was reached, " 

815 "so PICOS will ask CPLEX to stop the search.") 

816 self.abort() 

817 

818 # Register the callback handler with CPLEX. 

819 callback = self.int.register_callback(PicosInfoCallback) 

820 

821 # Pass parameters to the callback handler. Note that 

822 # callback.startTime will be set just before optimization begins. 

823 callback.printer = self._verbose 

824 callback.ubound = o.cplex_upr_bnd_limit 

825 callback.lbound = o.cplex_lwr_bnd_limit 

826 callback.bounds = [] if o.cplex_bnd_monitor else None 

827 

828 # Inform CPLEX about the problem type. 

829 # This seems necessary, as otherwise LP can get solved as MIP, producing 

830 # misleading status output (e.g. "not integer feasible"). 

831 conTypes = set(c.__class__ for c in self.ext.constraints.values()) 

832 quadObj = isinstance(self.ext.no.function, QuadraticExpression) 

833 cplexTypes = self.int.problem_type 

834 

835 if quadObj: 

836 if conTypes.issubset(set([DummyConstraint, AffineConstraint])): 

837 cplexType = cplexTypes.QP if continuous else cplexTypes.MIQP 

838 else: 

839 # Assume quadratic constraint types. 

840 cplexType = cplexTypes.QCP if continuous else cplexTypes.MIQCP 

841 else: 

842 if conTypes.issubset(set([DummyConstraint, AffineConstraint])): 

843 cplexType = cplexTypes.LP if continuous else cplexTypes.MILP 

844 else: 

845 # Assume quadratic constraint types. 

846 cplexType = cplexTypes.QCP if continuous else cplexTypes.MIQCP 

847 

848 # Silence a warning explaining that optimality target 3 changes the 

849 # problem type from QP to MIQP by doing so manually. 

850 if nonconvex_quad_obj: 

851 # Enforce consistency with CPLEXSolver.supports. 

852 assert cplexType in (cplexTypes.QP, cplexTypes.MIQP) 

853 

854 if p.optimalitytarget.get() == 3: # User might have changed it. 

855 cplexType = cplexTypes.MIQP 

856 

857 if cplexType is not None: 

858 self.int.set_problem_type(cplexType) 

859 

860 # Attempt to solve the problem. 

861 if callback: 

862 callback.startTime = time.time() 

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

864 try: 

865 if maxNumSolutions > 1: 

866 self.int.populate_solution_pool() 

867 numSolutions = self.int.solution.pool.get_num() 

868 else: 

869 self.int.solve() 

870 numSolutions = 1 

871 except cplex.exceptions.errors.CplexSolverError as error: 

872 if error.args[2] == 5002: 

873 self._handle_continuous_nonconvex_error(error) 

874 else: 

875 raise 

876 

877 solutions = [] 

878 for solutionNum in range(numSolutions): 

879 # Retrieve primals. 

880 primals = {} 

881 if o.primals is not False: 

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

883 try: 

884 indices = self._lookup(self._cplexVar, picosVar) 

885 

886 if maxNumSolutions > 1: 

887 value = self.int.solution.pool.get_values( 

888 solutionNum, indices) 

889 else: 

890 value = self.int.solution.get_values(indices) 

891 

892 primals[picosVar] = value 

893 except cplex.exceptions.errors.CplexSolverError: 

894 primals[picosVar] = None 

895 

896 # Retrieve duals. 

897 duals = {} 

898 if o.duals is not False and continuous: 

899 assert maxNumSolutions == 1 

900 

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

902 if isinstance(picosCon, DummyConstraint): 

903 duals[picosCon] = cvxopt.spmatrix( 

904 [], [], [], picosCon.size) 

905 continue 

906 

907 try: 

908 if isinstance(picosCon, AffineConstraint): 

909 indices = self._lookup(self._cplexLinCon, picosCon) 

910 values = self.int.solution.get_dual_values(indices) 

911 picosDual = cvxopt.matrix(values, picosCon.size) 

912 

913 if not picosCon.is_increasing(): 

914 picosDual = -picosDual 

915 elif isinstance(picosCon, SOCConstraint): 

916 metaCon = self._cplexMetaCon[picosCon] 

917 indices = self._lookup(self._cplexLinCon, metaCon) 

918 values = self.int.solution.get_dual_values(indices) 

919 picosDual = -cvxopt.matrix(values) 

920 picosDual[0] = -picosDual[0] 

921 elif isinstance(picosCon, RSOCConstraint): 

922 metaCon = self._cplexMetaCon[picosCon] 

923 indices = self._lookup(self._cplexLinCon, metaCon) 

924 values = self.int.solution.get_dual_values(indices) 

925 picosDual = -cvxopt.matrix(values) 

926 picosDual[0] = -picosDual[0] 

927 picosDual[1] = -picosDual[1] 

928 elif isinstance(picosCon, ConvexQuadraticConstraint): 

929 picosDual = None 

930 else: 

931 assert False, "Unexpected constraint type." 

932 

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

934 picosDual = -picosDual 

935 except cplex.exceptions.errors.CplexSolverError: 

936 duals[picosCon] = None 

937 else: 

938 duals[picosCon] = picosDual 

939 

940 # Retrieve objective value. 

941 try: 

942 if quadObj: 

943 # FIXME: Retrieval of QP and MIQP objective value appears to 

944 # miss the quadratic part. 

945 value = None 

946 elif maxNumSolutions > 1: 

947 value = self.int.solution.pool.get_objective_value( 

948 solutionNum) 

949 else: 

950 value = self.int.solution.get_objective_value() 

951 except cplex.exceptions.errors.CplexSolverError: 

952 value = None 

953 

954 # Retrieve solution status. 

955 code = self.int.solution.get_status() 

956 if code in CPLEX_STATUS_CODES: 

957 prmlStatus, dualStatus, probStatus = CPLEX_STATUS_CODES[code] 

958 else: 

959 prmlStatus = SS_UNKNOWN 

960 dualStatus = SS_UNKNOWN 

961 probStatus = PS_UNKNOWN 

962 

963 info = {} 

964 if o.cplex_bnd_monitor: 

965 info["bounds_monitor"] = callback.bounds 

966 

967 solutions.append(self._make_solution(value, primals, duals, 

968 prmlStatus, dualStatus, probStatus, info)) 

969 

970 if maxNumSolutions > 1: 

971 return solutions 

972 else: 

973 assert len(solutions) == 1 

974 return solutions[0] 

975 

976 

977# -------------------------------------- 

978__all__ = api_end(_API_START, globals())