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:`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 # TODO: Allow nonconvex quadratic constraints in the integer case? 

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

108 # the problem type. 

109 SUPPORTED = Specification( 

110 objectives=[ 

111 AffineExpression, 

112 QuadraticExpression], 

113 constraints=[ 

114 DummyConstraint, 

115 AffineConstraint, 

116 SOCConstraint, 

117 RSOCConstraint, 

118 ConvexQuadraticConstraint]) 

119 

120 @classmethod 

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

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

123 result = Solver.supports(footprint, explain) 

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

125 return result 

126 

127 if footprint.nonconvex_quadratic_objective and footprint.continuous: 

128 if explain: 

129 return (False, 

130 "Continuous problems with nonconvex quadratic objectives.") 

131 else: 

132 return False 

133 

134 if footprint not in cls.SUPPORTED: 

135 if explain: 

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

137 else: 

138 return False 

139 

140 return (True, None) if explain else True 

141 

142 @classmethod 

143 def default_penalty(cls): 

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

145 return 0.0 # Commercial solver. 

146 

147 @classmethod 

148 def test_availability(cls): 

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

150 cls.check_import("cplex") 

151 

152 @classmethod 

153 def names(cls): 

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

155 return "cplex", "CPLEX", "IBM ILOG CPLEX Optimization Studio" 

156 

157 @classmethod 

158 def is_free(cls): 

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

160 return False 

161 

162 CplexSOCC = namedtuple("CplexSOCC", 

163 ("LHSVars", "RHSVar", "LHSCons", "RHSCon", "quadCon")) 

164 

165 CplexRSOCC = namedtuple("CplexRSOCC", 

166 ("LHSVars", "RHSVars", "LHSCons", "RHSCons", "quadCon")) 

167 

168 def __init__(self, problem): 

169 """Initialize a CPLEX solver interface. 

170 

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

172 """ 

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

174 

175 self._cplexVarName = dict() 

176 """Maps PICOS variable indices to CPLEX variable names.""" 

177 

178 self._cplexLinConNames = dict() 

179 """Maps a PICOS (multidimensional) linear constraint to a collection of 

180 CPLEX (scalar) linear constraint names.""" 

181 

182 self._cplexQuadConName = dict() 

183 """Maps a PICOS quadratic or conic quadratic constraint to a CPLEX 

184 quadratic constraint name.""" 

185 

186 self._cplexSOCC = dict() 

187 """Maps a PICOS second order cone constraint to its CPLEX representation 

188 involving auxiliary variables and constraints.""" 

189 

190 self._cplexRSOCC = dict() 

191 """Maps a PICOS rotated second order cone constraint to its CPLEX 

192 representation involving auxiliary variables and constraints.""" 

193 

194 self.nextConstraintID = 0 

195 """Used to create unique names for constraints.""" 

196 

197 def __del__(self): 

198 if self.int is not None: 

199 self.int.end() 

200 

201 def reset_problem(self): 

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

203 if self.int is not None: 

204 self.int.end() 

205 self.int = None 

206 self._cplexVarName.clear() 

207 self._cplexLinConNames.clear() 

208 self._cplexQuadConName.clear() 

209 self._cplexSOCC.clear() 

210 self._cplexRSOCC.clear() 

211 

212 def _get_unique_constraint_id(self): 

213 ID = self.nextConstraintID 

214 self.nextConstraintID += 1 

215 return ID 

216 

217 def _make_cplex_var_names(self, picosVar, localIndex=None): 

218 """Make CPLEX variable names. 

219 

220 Converts a PICOS variable to a list of CPLEX variable names, each 

221 corresponding to one scalar variable contained in the PICOS variable. 

222 If localIndex is given, then only the name of the CPLEX variable 

223 representing the scalar variable with that offset is returned. 

224 The name format is "picosName[localIndex]". 

225 """ 

226 # TODO: This function appears in multiple solvers, move it to the Solver 

227 # base class as "_make_scalar_var_names". 

228 if localIndex is not None: 

229 return "{}[{}]".format(picosVar.name, localIndex) 

230 else: 

231 return [ 

232 self._make_cplex_var_names(picosVar, localIndex) 

233 for localIndex in range(picosVar.dim)] 

234 

235 def _import_variable(self, picosVar): 

236 import cplex 

237 

238 dim = picosVar.dim 

239 

240 # Create names. 

241 names = self._make_cplex_var_names(picosVar) 

242 

243 # Retrieve types. 

244 if isinstance(picosVar, CONTINUOUS_VARTYPES): 

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

246 elif isinstance(picosVar, IntegerVariable): 

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

248 elif isinstance(picosVar, BinaryVariable): 

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

250 else: 

251 assert False, "Unexpected variable type." 

252 

253 # Retrieve bounds. 

254 lowerBounds = [-cplex.infinity]*dim 

255 upperBounds = [cplex.infinity]*dim 

256 lower, upper = picosVar.bound_dicts 

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

258 lowerBounds[i] = b 

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

260 upperBounds[i] = b 

261 

262 # Import variable. 

263 # Note that CPLEX allows importing the objective function coefficients 

264 # for the new variables here, but that is done later to streamline 

265 # updates to the objective. 

266 self.int.variables.add( 

267 lb=lowerBounds, ub=upperBounds, types=types, names=names) 

268 

269 # Map PICOS indices to CPLEX names. 

270 for localIndex in range(dim): 

271 self._cplexVarName[picosVar.id_at(localIndex)] = names[localIndex] 

272 

273 if self._debug(): 

274 cplexVar = {"names": names, "types": types, 

275 "lowerBounds": lowerBounds, "upperBounds": upperBounds} 

276 self._debug( 

277 "Variable imported: {} → {}".format(picosVar, cplexVar)) 

278 

279 def _remove_variable(self, picosVar): 

280 cplexVarNames = [self._cplexVarName.pop(picosVar.id_at(localIndex)) 

281 for localIndex in range(picosVar.dim)] 

282 self.int.variables.delete(cplexVarNames) 

283 

284 def _affinexp_pic2cpl(self, picosExpression): 

285 import cplex 

286 for names, coefficients, constant in picosExpression.sparse_rows( 

287 None, indexFunction=self._make_cplex_var_names): 

288 yield cplex.SparsePair(ind=names, val=coefficients), constant 

289 

290 def _scalar_affinexp_pic2cpl(self, picosExpression): 

291 assert len(picosExpression) == 1 

292 return next(self._affinexp_pic2cpl(picosExpression)) 

293 

294 def _quadexp_pic2cpl(self, picosExpression): 

295 """Transform a quadratic expression from PICOS to CPLEX. 

296 

297 :returns: :class:`SparseTriple <cplex.SparseTriple>` mapping a pair of 

298 CPLEX variable names to scalar constants. 

299 """ 

300 import cplex 

301 

302 assert isinstance(picosExpression, QuadraticExpression) 

303 

304 cplexI, cplexJ, cplexV = [], [], [] 

305 for (picosVar1, picosVar2), picosCoefficients \ 

306 in picosExpression._sparse_quads.items(): 

307 for sparseIndex in range(len(picosCoefficients)): 

308 localVar1Index = picosCoefficients.I[sparseIndex] 

309 localVar2Index = picosCoefficients.J[sparseIndex] 

310 localCoefficient = picosCoefficients.V[sparseIndex] 

311 cplexI.append(self._cplexVarName[picosVar1.id + localVar1Index]) 

312 cplexJ.append(self._cplexVarName[picosVar2.id + localVar2Index]) 

313 cplexV.append(localCoefficient) 

314 

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

316 

317 def _import_linear_constraint(self, picosConstraint): 

318 import cplex 

319 

320 assert isinstance(picosConstraint, AffineConstraint) 

321 

322 length = len(picosConstraint) 

323 

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

325 cplexLHS, cplexRHS = [], [] 

326 for names, coefficients, constant in picosConstraint.sparse_Ab_rows( 

327 None, indexFunction=self._make_cplex_var_names): 

328 cplexLHS.append(cplex.SparsePair(ind=names, val=coefficients)) 

329 cplexRHS.append(constant) 

330 

331 # Retrieve senses. 

332 if picosConstraint.is_increasing(): 

333 senses = length * "L" 

334 elif picosConstraint.is_decreasing(): 

335 senses = length * "G" 

336 elif picosConstraint.is_equality(): 

337 senses = length * "E" 

338 else: 

339 assert False, "Unexpected constraint relation." 

340 

341 # Give unique names that are used to identify the constraint. This is 

342 # necessary as constraint indices can change if the problem is modified. 

343 conID = self._get_unique_constraint_id() 

344 names = ["{}:{}".format(conID, localConstraintIndex) 

345 for localConstraintIndex in range(length)] 

346 

347 if self._debug(): 

348 cplexConstraint = {"lin_expr": cplexLHS, "senses": senses, 

349 "rhs": cplexRHS, "names": names} 

350 self._debug( 

351 "Linear constraint imported: {} → {}".format( 

352 picosConstraint, cplexConstraint)) 

353 

354 # Import the constraint. 

355 self.int.linear_constraints.add( 

356 lin_expr=cplexLHS, senses=senses, rhs=cplexRHS, names=names) 

357 

358 return names 

359 

360 def _import_quad_constraint(self, picosConstraint): 

361 assert isinstance(picosConstraint, ConvexQuadraticConstraint) 

362 

363 # Retrieve 

364 # - CPLEX' linear term, which is the linear term of the affine 

365 # expression in the PICOS constraint, and 

366 # - CPLEX' right hand side, which is the negated constant term of the 

367 # affine expression in the PICOS constraint. 

368 cplexLinear, cplexRHS = \ 

369 self._scalar_affinexp_pic2cpl(picosConstraint.le0.aff) 

370 cplexRHS = -cplexRHS 

371 

372 # Retrieve CPLEX' quadratic term. 

373 cplexQuad = self._quadexp_pic2cpl(picosConstraint.le0) 

374 

375 # Give a unique name that is used to identify the constraint. This is 

376 # necessary as constraint indices can change if the problem is modified. 

377 name = "{}:{}".format(self._get_unique_constraint_id(), 0) 

378 

379 if self._debug(): 

380 cplexConstraint = {"lin_expr": cplexLinear, "quad_expr": cplexQuad, 

381 "rhs": cplexRHS, "name": name} 

382 self._debug( 

383 "Quadratic constraint imported: {} → {}".format( 

384 picosConstraint, cplexConstraint)) 

385 

386 # Import the constraint. 

387 self.int.quadratic_constraints.add(lin_expr=cplexLinear, 

388 quad_expr=cplexQuad, sense="L", rhs=cplexRHS, name=name) 

389 

390 return name 

391 

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

393 def _import_socone_constraint(self, picosConstraint): 

394 import cplex 

395 

396 assert isinstance(picosConstraint, SOCConstraint) 

397 

398 picosLHS = picosConstraint.ne 

399 picosRHS = picosConstraint.ub 

400 picosLHSLen = len(picosLHS) 

401 

402 # Make identifying names for the auxiliary variables and constraints. 

403 conID = self._get_unique_constraint_id() 

404 cplexLHSVars = ["{}:V{}".format(conID, i) for i in range(picosLHSLen)] 

405 cplexRHSVar = "{}:V{}".format(conID, picosLHSLen) 

406 cplexLHSCons = ["{}:C{}".format(conID, i) for i in range(picosLHSLen)] 

407 cplexRHSCon = "{}:C{}".format(conID, picosLHSLen) 

408 cplexQuadCon = "{}:C{}".format(conID, picosLHSLen + 1) 

409 

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

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

412 self.int.variables.add( 

413 names=cplexLHSVars, lb=[-cplex.infinity] * picosLHSLen, 

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

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

416 self.int.variables.add( 

417 names=[cplexRHSVar], lb=[0.0], ub=[+cplex.infinity], 

418 types=self.int.variables.type.continuous) 

419 

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

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

422 cplexLHSConsLHSs = [] 

423 cplexLHSConsRHSs = [] 

424 for localConIndex, (localLinExp, localConstant) in \ 

425 enumerate(self._affinexp_pic2cpl(picosLHS)): 

426 localConstant = -localConstant 

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

428 localLinExp.val.append(-1.0) 

429 cplexLHSConsLHSs.append(localLinExp) 

430 cplexLHSConsRHSs.append(localConstant) 

431 self.int.linear_constraints.add( 

432 names=cplexLHSCons, lin_expr=cplexLHSConsLHSs, 

433 senses="E" * picosLHSLen, rhs=cplexLHSConsRHSs) 

434 

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

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

437 cplexRHSConLHS, cplexRHSConRHS = \ 

438 self._scalar_affinexp_pic2cpl(-picosRHS) 

439 cplexRHSConRHS = -cplexRHSConRHS 

440 cplexRHSConLHS.ind.append(cplexRHSVar) 

441 cplexRHSConLHS.val.append(1.0) 

442 self.int.linear_constraints.add( 

443 names=[cplexRHSCon], lin_expr=[cplexRHSConLHS], 

444 senses="E", rhs=[cplexRHSConRHS]) 

445 

446 # Add a quadratic constraint over the auxiliary variables that 

447 # represents the PICOS second order cone constraint itself. 

448 quadIndices = [cplexRHSVar] + list(cplexLHSVars) 

449 quadExpr = cplex.SparseTriple( 

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

451 self.int.quadratic_constraints.add( 

452 name=cplexQuadCon, quad_expr=quadExpr, sense="L", rhs=0.0) 

453 

454 cplexMetaCon = self.CplexSOCC(LHSVars=cplexLHSVars, RHSVar=cplexRHSVar, 

455 LHSCons=cplexLHSCons, RHSCon=cplexRHSCon, quadCon=cplexQuadCon) 

456 

457 if self._debug(): 

458 cplexCons = { 

459 "LHSs of LHS auxiliary equalities": cplexLHSConsLHSs, 

460 "RHSs of LHS auxiliary equalities": cplexLHSConsRHSs, 

461 "LHS of RHS auxiliary equality": cplexRHSConLHS, 

462 "RHS of RHS auxiliary equality": cplexRHSConRHS, 

463 "Non-positive quadratic term": quadExpr} 

464 self._debug( 

465 "SOcone constraint imported: {} → {}, {}".format( 

466 picosConstraint, cplexMetaCon, cplexCons)) 

467 

468 return cplexMetaCon 

469 

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

471 def _import_rscone_constraint(self, picosConstraint): 

472 import cplex 

473 

474 assert isinstance(picosConstraint, RSOCConstraint) 

475 

476 picosLHS = picosConstraint.ne 

477 picosRHS1 = picosConstraint.ub1 

478 picosRHS2 = picosConstraint.ub2 

479 picosLHSLen = len(picosLHS) 

480 

481 # Make identifying names for the auxiliary variables and constraints. 

482 conID = self._get_unique_constraint_id() 

483 cplexLHSVars = ["{}:V{}".format(conID, i) for i in range(picosLHSLen)] 

484 cplexRHSVars = ["{}:V{}".format(conID, picosLHSLen + i) for i in (0, 1)] 

485 cplexLHSCons = ["{}:C{}".format(conID, i) for i in range(picosLHSLen)] 

486 cplexRHSCons = ["{}:C{}".format(conID, picosLHSLen + i) for i in (0, 1)] 

487 cplexQuadCon = "{}:C{}".format(conID, picosLHSLen + 2) 

488 

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

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

491 self.int.variables.add( 

492 names=cplexLHSVars, lb=[-cplex.infinity] * picosLHSLen, 

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

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

495 self.int.variables.add( 

496 names=cplexRHSVars, lb=[0.0, 0.0], ub=[+cplex.infinity] * 2, 

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

498 

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

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

501 cplexLHSConsLHSs = [] 

502 cplexLHSConsRHSs = [] 

503 for localConIndex, (localLinExp, localConstant) in \ 

504 enumerate(self._affinexp_pic2cpl(picosLHS)): 

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

506 localLinExp.val.append(-1.0) 

507 localConstant = -localConstant 

508 cplexLHSConsLHSs.append(localLinExp) 

509 cplexLHSConsRHSs.append(localConstant) 

510 self.int.linear_constraints.add( 

511 names=cplexLHSCons, lin_expr=cplexLHSConsLHSs, 

512 senses="E" * picosLHSLen, rhs=cplexLHSConsRHSs) 

513 

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

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

516 cplexRHSConsLHSs = [] 

517 cplexRHSConsRHSs = [] 

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

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

520 linExp.ind.append(cplexRHSVar) 

521 linExp.val.append(1.0) 

522 constant = -constant 

523 cplexRHSConsLHSs.append(linExp) 

524 cplexRHSConsRHSs.append(constant) 

525 self.int.linear_constraints.add( 

526 names=cplexRHSCons, lin_expr=cplexRHSConsLHSs, 

527 senses="E" * 2, rhs=cplexRHSConsRHSs) 

528 

529 # Add a quadratic constraint over the auxiliary variables that 

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

531 quadExpr = cplex.SparseTriple( 

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

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

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

535 self.int.quadratic_constraints.add( 

536 name=cplexQuadCon, quad_expr=quadExpr, sense="L", rhs=0.0) 

537 

538 cplexMetaCon = self.CplexRSOCC( 

539 LHSVars=cplexLHSVars, RHSVars=cplexRHSVars, LHSCons=cplexLHSCons, 

540 RHSCons=cplexRHSCons, quadCon=cplexQuadCon) 

541 

542 if self._debug(): 

543 cplexCons = { 

544 "LHSs of LHS auxiliary equalities": cplexLHSConsLHSs, 

545 "RHSs of LHS auxiliary equalities": cplexLHSConsRHSs, 

546 "LHSs of RHS auxiliary equalities": cplexRHSConsLHSs, 

547 "RHSs of RHS auxiliary equalities": cplexRHSConsRHSs, 

548 "Non-positive quadratic term": quadExpr} 

549 self._debug( 

550 "RScone constraint imported: {} → {}, {}".format( 

551 picosConstraint, cplexMetaCon, cplexCons)) 

552 

553 return cplexMetaCon 

554 

555 def _import_constraint(self, picosConstraint): 

556 # Import constraint based on type and keep track of the corresponding 

557 # CPLEX constraint and auxiliary variable names. 

558 if isinstance(picosConstraint, AffineConstraint): 

559 self._cplexLinConNames[picosConstraint] = \ 

560 self._import_linear_constraint(picosConstraint) 

561 elif isinstance(picosConstraint, ConvexQuadraticConstraint): 

562 self._cplexQuadConName[picosConstraint] = \ 

563 self._import_quad_constraint(picosConstraint) 

564 elif isinstance(picosConstraint, SOCConstraint): 

565 self._cplexSOCC[picosConstraint] = \ 

566 self._import_socone_constraint(picosConstraint) 

567 elif isinstance(picosConstraint, RSOCConstraint): 

568 self._cplexRSOCC[picosConstraint] = \ 

569 self._import_rscone_constraint(picosConstraint) 

570 else: 

571 assert isinstance(picosConstraint, DummyConstraint), \ 

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

573 picosConstraint.__class__.__name__) 

574 

575 def _remove_constraint(self, picosConstraint): 

576 if isinstance(picosConstraint, AffineConstraint): 

577 self.int.linear_constraints.delete( 

578 self._cplexLinConNames.pop(picosConstraint)) 

579 elif isinstance(picosConstraint, ConvexQuadraticConstraint): 

580 self.int.quadratic_constraints.delete( 

581 self._cplexQuadConName.pop(picosConstraint)) 

582 elif isinstance(picosConstraint, SOCConstraint): 

583 c = self._cplexSOCC.pop(picosConstraint) 

584 self.int.linear_constraints.delete(c.cplexLHSCons + [c.cplexRHSCon]) 

585 self.int.quadratic_constraints.delete(c.cplexQuadCon) 

586 self.int.variables.delete(c.cplexLHSVars + [c.cplexRHSVar]) 

587 elif isinstance(picosConstraint, RSOCConstraint): 

588 c = self._cplexRSOCC.pop(picosConstraint) 

589 self.int.linear_constraints.delete(c.cplexLHSCons + c.cplexRHSCons) 

590 self.int.quadratic_constraints.delete(c.cplexQuadCon) 

591 self.int.variables.delete(c.cplexLHSVars + c.cplexRHSVars) 

592 else: 

593 assert isinstance(picosConstraint, DummyConstraint), \ 

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

595 picosConstraint.__class__.__name__) 

596 

597 def _import_affine_objective(self, picosExpression): 

598 assert isinstance(picosExpression, AffineExpression) 

599 

600 if picosExpression._constant_coef: 

601 offset = picosExpression._constant_coef[0] 

602 

603 self.int.objective.set_offset(offset) 

604 

605 if self._debug(): 

606 self._debug("Constant part of objective imported: {} → {}" 

607 .format(picosExpression.cst.string, offset)) 

608 

609 cplexExpression = [] 

610 for picosVar, picosCoefficient in picosExpression._linear_coefs.items(): 

611 assert picosCoefficient.size[0] == 1 

612 

613 for localIndex in range(picosVar.dim): 

614 cplexCoefficient = picosCoefficient[localIndex] 

615 if not cplexCoefficient: 

616 continue 

617 picosIndex = picosVar.id + localIndex 

618 cplexName = self._cplexVarName[picosIndex] 

619 cplexExpression.append((cplexName, cplexCoefficient)) 

620 

621 if cplexExpression: 

622 self.int.objective.set_linear(cplexExpression) 

623 

624 if self._debug(): 

625 self._debug("Linear part of objective imported: {} → {}" 

626 .format(picosExpression.lin.string, cplexExpression)) 

627 

628 def _reset_affine_objective(self): 

629 self.int.objective.set_offset(0.0) 

630 

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

632 if any(linear): 

633 self.int.objective.set_linear( 

634 [(cplexVarIndex, 0.0) for cplexVarIndex, coefficient 

635 in enumerate(linear) if coefficient]) 

636 

637 def _import_quadratic_objective(self, picosExpression): 

638 assert isinstance(picosExpression, QuadraticExpression) 

639 

640 # Import affine part of objective function. 

641 self._import_affine_objective(picosExpression.aff) 

642 

643 # Import quadratic part of objective function. 

644 cplexQuadExpression = self._quadexp_pic2cpl(picosExpression) 

645 cplexQuadCoefficients = zip( 

646 cplexQuadExpression.ind1, cplexQuadExpression.ind2, 

647 [2.0 * coefficient for coefficient in cplexQuadExpression.val]) 

648 self.int.objective.set_quadratic_coefficients(cplexQuadCoefficients) 

649 

650 self._debug("Quadratic part of objective imported: {} → {}" 

651 .format(picosExpression.quad.string, cplexQuadCoefficients)) 

652 

653 def _reset_quadratic_objective(self): 

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

655 if quadratics: 

656 self.int.objective.set_quadratic( 

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

658 for sparsePair in quadratics]) 

659 

660 def _import_objective(self): 

661 picosSense, picosObjective = self.ext.no 

662 

663 # Import objective sense. 

664 if picosSense == "min": 

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

666 else: 

667 assert picosSense == "max" 

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

669 self.int.objective.set_sense(cplexSense) 

670 

671 # Import objective function. 

672 if isinstance(picosObjective, AffineExpression): 

673 self._import_affine_objective(picosObjective) 

674 else: 

675 assert isinstance(picosObjective, QuadraticExpression) 

676 self._import_quadratic_objective(picosObjective) 

677 

678 def _reset_objective(self): 

679 self._reset_affine_objective() 

680 self._reset_quadratic_objective() 

681 

682 def _import_problem(self): 

683 import cplex 

684 

685 # Create a problem instance. 

686 self.int = cplex.Cplex() 

687 

688 # Import variables. 

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

690 self._import_variable(variable) 

691 

692 # Import constraints. 

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

694 self._import_constraint(constraint) 

695 

696 # Set objective. 

697 self._import_objective() 

698 

699 def _update_problem(self): 

700 for oldConstraint in self._removed_constraints(): 

701 self._remove_constraint(oldConstraint) 

702 

703 for oldVariable in self._removed_variables(): 

704 self._remove_variable(oldVariable) 

705 

706 for newVariable in self._new_variables(): 

707 self._import_variable(newVariable) 

708 

709 for newConstraint in self._new_constraints(): 

710 self._import_constraint(newConstraint) 

711 

712 if self._objective_has_changed(): 

713 self._reset_objective() 

714 self._import_objective() 

715 

716 def _solve(self): 

717 import cplex 

718 

719 # Reset options. 

720 self.int.parameters.reset() 

721 

722 o = self.ext.options 

723 p = self.int.parameters 

724 

725 continuous = self.ext.is_continuous() 

726 

727 # verbosity 

728 verbosity = self.verbosity() 

729 if verbosity <= 0: 

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

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

732 # option that is set. 

733 self.int.set_results_stream(None) 

734 else: 

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

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

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

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

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

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

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

742 

743 # abs_prim_fsb_tol 

744 if o.abs_prim_fsb_tol is not None: 

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

746 

747 # abs_dual_fsb_tol 

748 if o.abs_dual_fsb_tol is not None: 

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

750 

751 # rel_prim_fsb_tol, rel_dual_fsb_tol, rel_ipm_opt_tol 

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

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

754 if convergenceTols: 

755 convergenceTol = min(convergenceTols) 

756 p.barrier.convergetol.set(convergenceTol) 

757 p.barrier.qcpconvergetol.set(convergenceTol) 

758 

759 # abs_bnb_opt_tol 

760 if o.abs_bnb_opt_tol is not None: 

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

762 

763 # rel_bnb_opt_tol 

764 if o.rel_bnb_opt_tol is not None: 

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

766 

767 # integrality_tol 

768 if o.integrality_tol is not None: 

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

770 

771 # markowitz_tol 

772 if o.markowitz_tol is not None: 

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

774 

775 # max_iterations 

776 if o.max_iterations is not None: 

777 maxit = o.max_iterations 

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

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

780 

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

782 

783 # lp_node_method 

784 if o.lp_node_method is not None: 

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

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

787 

788 # lp_root_method 

789 if o.lp_root_method is not None: 

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

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

792 

793 # timelimit 

794 if o.timelimit is not None: 

795 p.timelimit.set(o.timelimit) 

796 

797 # treememory 

798 if o.treememory is not None: 

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

800 

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

802 if o.max_fsb_nodes is not None \ 

803 and o.pool_size is not None: 

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

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

806 

807 # max_fsb_nodes 

808 if o.max_fsb_nodes is not None: 

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

810 

811 # pool_size 

812 if o.pool_size is not None: 

813 if continuous: 

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

815 "be used with mixed integer problems.") 

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

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

818 else: 

819 maxNumSolutions = 1 

820 

821 # pool_relgap 

822 if o.pool_rel_gap is not None: 

823 if o.pool_size is None: 

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

825 "the option 'pool_size'.") 

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

827 

828 # pool_abs_gap 

829 if o.pool_abs_gap is not None: 

830 if o.pool_size is None: 

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

832 "the option 'pool_size'.") 

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

834 

835 # hotstart 

836 if o.hotstart: 

837 names, values = [], [] 

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

839 if picosVar.valued: 

840 for localIndex in range(picosVar.dim): 

841 name = self._cplexVarName[picosVar.id_at(localIndex)] 

842 names.append(name) 

843 values.append(picosVar.internal_value[localIndex]) 

844 if names: 

845 self.int.MIP_starts.add( 

846 cplex.SparsePair(ind=names, val=values), 

847 self.int.MIP_starts.effort_level.repair) 

848 

849 # Handle CPLEX-specific options. 

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

851 try: 

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

853 except AttributeError as error: 

854 self._handle_bad_solver_specific_option_key(key, error) 

855 

856 try: 

857 parameter.set(value) 

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

859 self._handle_bad_solver_specific_option_value(key, value, error) 

860 

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

862 # "cplex_bnd_monitor" via a CPLEX callback handler. 

863 callback = None 

864 if o.cplex_upr_bnd_limit or o.cplex_lwr_bnd_limit \ 

865 or o.cplex_bnd_monitor: 

866 from cplex.callbacks import MIPInfoCallback 

867 

868 class PicosInfoCallback(MIPInfoCallback): 

869 def __call__(self): 

870 v1 = self.get_incumbent_objective_value() 

871 v2 = self.get_best_objective_value() 

872 ub = max(v1, v2) 

873 lb = min(v1, v2) 

874 if self.bounds is not None: 

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

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

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

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

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

880 self.abort() 

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

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

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

884 self.abort() 

885 

886 # Register the callback handler with CPLEX. 

887 callback = self.int.register_callback(PicosInfoCallback) 

888 

889 # Pass parameters to the callback handler. Note that 

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

891 callback.printer = self._verbose 

892 callback.ubound = o.cplex_upr_bnd_limit 

893 callback.lbound = o.cplex_lwr_bnd_limit 

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

895 

896 # Inform CPLEX about the problem type. 

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

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

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

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

901 cplexTypes = self.int.problem_type 

902 

903 if quadObj: 

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

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

906 else: 

907 # Assume quadratic constraint types. 

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

909 else: 

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

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

912 else: 

913 # Assume quadratic constraint types. 

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

915 

916 if cplexType is not None: 

917 self.int.set_problem_type(cplexType) 

918 

919 # Attempt to solve the problem. 

920 if callback: 

921 callback.startTime = time.time() 

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

923 try: 

924 if maxNumSolutions > 1: 

925 self.int.populate_solution_pool() 

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

927 else: 

928 self.int.solve() 

929 numSolutions = 1 

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

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

932 self._handle_continuous_nonconvex_error(error) 

933 else: 

934 raise 

935 

936 solutions = [] 

937 for solutionNum in range(numSolutions): 

938 # Retrieve primals. 

939 primals = {} 

940 if o.primals is not False: 

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

942 try: 

943 cplexNames = [] 

944 

945 for localIndex in range(picosVar.dim): 

946 picosIndex = picosVar.id + localIndex 

947 cplexNames.append(self._cplexVarName[picosIndex]) 

948 if maxNumSolutions > 1: 

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

950 solutionNum, cplexNames) 

951 else: 

952 value = self.int.solution.get_values(cplexNames) 

953 primals[picosVar] = value 

954 except cplex.exceptions.errors.CplexSolverError: 

955 primals[picosVar] = None 

956 

957 # Retrieve duals. 

958 duals = {} 

959 if o.duals is not False and continuous: 

960 assert maxNumSolutions == 1 

961 

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

963 if isinstance(picosCon, DummyConstraint): 

964 duals[picosCon] = cvxopt.spmatrix( 

965 [], [], [], picosCon.size) 

966 continue 

967 

968 try: 

969 if isinstance(picosCon, AffineConstraint): 

970 cplexCons = self._cplexLinConNames[picosCon] 

971 values = self.int.solution.get_dual_values( 

972 cplexCons) 

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

974 

975 # Flip sign based on constraint relation. 

976 if not picosCon.is_increasing(): 

977 picosDual = -picosDual 

978 elif isinstance(picosCon, SOCConstraint): 

979 cplexMetaCon = self._cplexSOCC[picosCon] 

980 lb = self.int.solution.get_dual_values( 

981 cplexMetaCon.RHSCon) 

982 z = self.int.solution.get_dual_values( 

983 list(cplexMetaCon.LHSCons)) 

984 picosDual = -cvxopt.matrix([-lb] + z) 

985 elif isinstance(picosCon, RSOCConstraint): 

986 cplexMetaCon = self._cplexRSOCC[picosCon] 

987 ab = [-x for x in self.int.solution.get_dual_values( 

988 list(cplexMetaCon.RHSCons))] 

989 z = self.int.solution.get_dual_values( 

990 list(cplexMetaCon.LHSCons)) 

991 picosDual = -cvxopt.matrix(ab + z) 

992 elif isinstance(picosCon, ConvexQuadraticConstraint): 

993 picosDual = None 

994 else: 

995 assert False, "Unexpected constraint type." 

996 

997 # Flip sign based on objective sense. 

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

999 picosDual = -picosDual 

1000 except cplex.exceptions.errors.CplexSolverError: 

1001 duals[picosCon] = None 

1002 else: 

1003 duals[picosCon] = picosDual 

1004 

1005 # Retrieve objective value. 

1006 try: 

1007 if quadObj: 

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

1009 # miss the quadratic part. 

1010 value = None 

1011 elif maxNumSolutions > 1: 

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

1013 solutionNum) 

1014 else: 

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

1016 except cplex.exceptions.errors.CplexSolverError: 

1017 value = None 

1018 

1019 # Retrieve solution status. 

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

1021 if code in CPLEX_STATUS_CODES: 

1022 prmlStatus, dualStatus, probStatus = CPLEX_STATUS_CODES[code] 

1023 else: 

1024 prmlStatus = SS_UNKNOWN 

1025 dualStatus = SS_UNKNOWN 

1026 probStatus = PS_UNKNOWN 

1027 

1028 info = {} 

1029 if o.cplex_bnd_monitor: 

1030 info["bounds_monitor"] = callback.bounds 

1031 

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

1033 prmlStatus, dualStatus, probStatus, info)) 

1034 

1035 if maxNumSolutions > 1: 

1036 return solutions 

1037 else: 

1038 assert len(solutions) == 1 

1039 return solutions[0] 

1040 

1041 

1042# -------------------------------------- 

1043__all__ = api_end(_API_START, globals())