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

20 

21import cvxopt 

22 

23from ..apidoc import api_end, api_start 

24from ..constraints import (AffineConstraint, ConvexQuadraticConstraint, 

25 DummyConstraint, NonconvexQuadraticConstraint, 

26 RSOCConstraint, SOCConstraint) 

27from ..expressions import AffineExpression, BinaryVariable, IntegerVariable 

28from ..modeling.footprint import Specification 

29from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED, 

30 PS_UNKNOWN, SS_INFEASIBLE, SS_OPTIMAL, 

31 SS_PREMATURE, SS_UNKNOWN) 

32from .solver import OptionValueError, Solver 

33 

34_API_START = api_start(globals()) 

35# ------------------------------- 

36 

37 

38class SCIPSolver(Solver): 

39 """Interface to the SCIP solver via the PySCIPOpt Python interface.""" 

40 

41 SUPPORTED = Specification( 

42 objectives=[ 

43 AffineExpression], 

44 constraints=[ 

45 DummyConstraint, 

46 AffineConstraint, 

47 SOCConstraint, 

48 RSOCConstraint, 

49 ConvexQuadraticConstraint, 

50 NonconvexQuadraticConstraint]) 

51 

52 @classmethod 

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

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

55 result = Solver.supports(footprint, explain) 

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

57 return result 

58 

59 # HACK: SCIP has trouble with continuous problems, in particular it only 

60 # produces dual values for LPs (wehich are sometimes incorrect). 

61 # Nonconvex quadratics are supported as no other solver can deal 

62 # with them at this point. 

63 if footprint.continuous \ 

64 and not ("con", NonconvexQuadraticConstraint) in footprint: 

65 if explain: 

66 return (False, 

67 "Continuous problems, except for nonconvex quadratic ones " 

68 "(PICOS' choice).") 

69 else: 

70 return False 

71 

72 if footprint not in cls.SUPPORTED: 

73 if explain: 

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

75 else: 

76 return False 

77 

78 return (True, None) if explain else True 

79 

80 @classmethod 

81 def default_penalty(cls): 

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

83 return 1.0 # Stable free/open source solver. 

84 

85 @classmethod 

86 def test_availability(cls): 

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

88 cls.check_import("pyscipopt") 

89 

90 @classmethod 

91 def names(cls): 

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

93 return "scip", "SCIP", "SCIP Optimization Suite" 

94 

95 @classmethod 

96 def is_free(cls): 

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

98 return False 

99 

100 def __init__(self, problem): 

101 """Initialize a SCIP solver interface. 

102 

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

104 """ 

105 super(SCIPSolver, self).__init__(problem) 

106 

107 self._scipVar = dict() 

108 """Maps PICOS variable indices to SCIP variables.""" 

109 

110 self._scipCons = dict() 

111 """ 

112 Maps PICOS constraints to lists of SCIP constraints. 

113 

114 For PICOS second order cone constraints, the first entry in the list is 

115 a SCIP quadratic constraint and the second entry is a SCIP linear 

116 auxiliary constraint. 

117 """ 

118 

119 self._scipQuadObjAuxVar = None 

120 """A SCIP auxiliary variable to support a PICOS quadratic objective.""" 

121 

122 self._scipQuadObjAuxCon = None 

123 """A SCIP auxiliary constraint to support a PICOS quadr. objective.""" 

124 

125 def reset_problem(self): 

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

127 self.int = None 

128 

129 self._scipVar.clear() 

130 self._scipCons.clear() 

131 

132 self._scipQuadObjAuxVar = None 

133 self._scipQuadObjAuxCon = None 

134 

135 def _make_scip_var_names(self, picosVar, localIndex=None): 

136 """Make SCIP variable names. 

137 

138 Converts a PICOS variable to a list of SCIP variable names, each 

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

140 If localIndex is given, then only the name of the SCIP variable 

141 representing the scalar variable with that offset is returned. 

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

143 """ 

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

145 # base class as "_make_scalar_var_names". 

146 if localIndex is not None: 

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

148 else: 

149 return [ 

150 self._make_scip_var_names(picosVar, localIndex) 

151 for localIndex in range(len(picosVar))] 

152 

153 def _import_variable(self, picosVar): 

154 from math import ceil, floor 

155 

156 dim = picosVar.dim 

157 inf = self.int.infinity() 

158 

159 # Retrieve type code. 

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

161 scipVarType = "I" 

162 else: 

163 scipVarType = "C" 

164 

165 # Retrieve bounds. 

166 lowerBounds = [-inf]*dim 

167 upperBounds = [inf]*dim 

168 lower, upper = picosVar.bound_dicts 

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

170 lowerBounds[i] = b 

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

172 upperBounds[i] = b 

173 

174 # Refine bounds. 

175 if isinstance(picosVar, IntegerVariable): 

176 lowerBounds = [int(ceil(b)) for b in lowerBounds] 

177 upperBounds = [int(floor(b)) for b in upperBounds] 

178 

179 # Give names. 

180 scipNames = self._make_scip_var_names(picosVar) 

181 

182 # Import scalar variables. 

183 for i in range(dim): 

184 self._scipVar[picosVar.id_at(i)] = self.int.addVar( 

185 name=scipNames[i], vtype=scipVarType, 

186 lb=lowerBounds[i], ub=upperBounds[i]) 

187 

188 def _import_variable_values(self): 

189 # TODO: Import variable values with SCIP. 

190 raise NotImplementedError 

191 

192 def _reset_variable_values(self): 

193 # TODO: Import variable values with SCIP. 

194 raise NotImplementedError 

195 

196 def _affinexp_pic2scip(self, picosExpression): 

197 """Transform an affine expression from PICOS to SCIP. 

198 

199 Multidimensional PICOS expressions are converted to a number of scalar 

200 SCIP expressions. 

201 

202 :returns: A :class:`list` of :class:`SCIP expressions 

203 <pyscipopt.scip.Expr>`. 

204 """ 

205 from pyscipopt.scip import Expr 

206 

207 if picosExpression is None: 

208 yield Expr() 

209 return 

210 

211 for scipVars, coefficients, constant in picosExpression.sparse_rows( 

212 None, indexFunction=lambda picosVar, localIndex: 

213 self._scipVar[picosVar.id_at(localIndex)]): 

214 scipExpression = Expr() 

215 for scipVar, coefficient in zip(scipVars, coefficients): 

216 scipExpression += coefficient * scipVar 

217 scipExpression += constant 

218 yield scipExpression 

219 

220 def _scalar_affinexp_pic2scip(self, picosExpression): 

221 """Transform a PICOS scalar affine expression into a SCIP expression. 

222 

223 :returns: A :class:`SCIP expression <pyscipopt.scip.Expr>`. 

224 """ 

225 assert len(picosExpression) == 1 

226 return next(self._affinexp_pic2scip(picosExpression)) 

227 

228 def _quadexp_pic2scip(self, picosExpression): 

229 # Convert affine part. 

230 scipExpression = self._scalar_affinexp_pic2scip(picosExpression.aff) 

231 

232 # Convert quadratic part. 

233 for (pVar1, pVar2), pCoefficients \ 

234 in picosExpression._sparse_quads.items(): 

235 for sparseIndex in range(len(pCoefficients)): 

236 localVar1Index = pCoefficients.I[sparseIndex] 

237 localVar2Index = pCoefficients.J[sparseIndex] 

238 localCoefficient = pCoefficients.V[sparseIndex] 

239 scipVar1 = self._scipVar[pVar1.id_at(localVar1Index)] 

240 scipVar2 = self._scipVar[pVar2.id_at(localVar2Index)] 

241 scipExpression += localCoefficient * scipVar1 * scipVar2 

242 

243 return scipExpression 

244 

245 def _import_linear_constraint(self, picosConstraint): 

246 assert isinstance(picosConstraint, AffineConstraint) 

247 

248 scipCons = [] 

249 picosLHS = picosConstraint.lhs - picosConstraint.rhs 

250 for scipLHS in self._affinexp_pic2scip(picosLHS): 

251 if picosConstraint.is_increasing(): 

252 scipCons.append(self.int.addCons(scipLHS <= 0.0)) 

253 elif picosConstraint.is_decreasing(): 

254 scipCons.append(self.int.addCons(scipLHS >= 0.0)) 

255 elif picosConstraint.is_equality(): 

256 scipCons.append(self.int.addCons(scipLHS == 0.0)) 

257 else: 

258 assert False, "Unexpected constraint relation." 

259 

260 return scipCons 

261 

262 def _import_quad_constraint(self, picosConstraint): 

263 assert isinstance(picosConstraint, 

264 (ConvexQuadraticConstraint, NonconvexQuadraticConstraint)) 

265 

266 scipLHS = self._quadexp_pic2scip(picosConstraint.le0) 

267 scipCons = [self.int.addCons(scipLHS <= 0.0)] 

268 

269 return scipCons 

270 

271 def _import_socone_constraint(self, picosConstraint): 

272 scipCons = [] 

273 scipQuadLHS = self._quadexp_pic2scip( 

274 (picosConstraint.ne | picosConstraint.ne) 

275 - (picosConstraint.ub * picosConstraint.ub)) 

276 scipCons.append(self.int.addCons(scipQuadLHS <= 0.0)) 

277 

278 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub) 

279 if scipAuxLHS.degree() > 0: 

280 scipCons.append(self.int.addCons(scipAuxLHS >= 0.0)) 

281 

282 return scipCons 

283 

284 def _import_rscone_constraint(self, picosConstraint): 

285 scipCons = [] 

286 scipLHS = self._quadexp_pic2scip( 

287 (picosConstraint.ne | picosConstraint.ne) 

288 - (picosConstraint.ub1 * picosConstraint.ub2)) 

289 scipCons.append(self.int.addCons(scipLHS <= 0.0)) 

290 

291 # Make sure that either the RHS is constant, or one of the two 

292 # expressions is non-negative. 

293 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub1) 

294 if scipAuxLHS.degree() > 0: 

295 scipCons.append(self.int.addCons(scipAuxLHS >= 0.0)) 

296 else: 

297 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub2) 

298 if scipAuxLHS.degree() > 0: 

299 scipCons.append(self.int.addCons(scipAuxLHS >= 0.0)) 

300 

301 return scipCons 

302 

303 def _import_constraint(self, picosConstraint): 

304 # Import constraint based on type. 

305 if isinstance(picosConstraint, AffineConstraint): 

306 scipCons = self._import_linear_constraint(picosConstraint) 

307 elif isinstance(picosConstraint, 

308 (ConvexQuadraticConstraint, NonconvexQuadraticConstraint)): 

309 scipCons = self._import_quad_constraint(picosConstraint) 

310 elif isinstance(picosConstraint, SOCConstraint): 

311 scipCons = self._import_socone_constraint(picosConstraint) 

312 elif isinstance(picosConstraint, RSOCConstraint): 

313 scipCons = self._import_rscone_constraint(picosConstraint) 

314 else: 

315 assert isinstance(picosConstraint, DummyConstraint), \ 

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

317 picosConstraint.__class__.__name__) 

318 

319 # Map PICOS constraints to lists of SCIP constraints. 

320 self._scipCons[picosConstraint] = scipCons 

321 

322 def _import_objective(self): 

323 picosSense, picosObjective = self.ext.no 

324 

325 # Retrieve objective sense. 

326 if picosSense == "min": 

327 scipSense = "minimize" 

328 else: 

329 assert picosSense == "max" 

330 scipSense = "maximize" 

331 

332 # Import objective function. 

333 scipObjective = self._scalar_affinexp_pic2scip(picosObjective) 

334 

335 self.int.setObjective(scipObjective, scipSense) 

336 

337 def _import_problem(self): 

338 import pyscipopt as scip 

339 

340 # Create a problem instance. 

341 self.int = scip.Model() 

342 

343 # Import variables. 

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

345 self._import_variable(variable) 

346 

347 # Import constraints. 

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

349 self._import_constraint(constraint) 

350 

351 # Set objective. 

352 self._import_objective() 

353 

354 def _update_problem(self): 

355 # TODO: Support all problem updates supported by SCIP. 

356 raise NotImplementedError 

357 

358 def _solve(self): 

359 import pyscipopt as scip 

360 

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

362 isLP = isinstance(self.ext.no.function, AffineExpression) \ 

363 and all(isinstance(constraint, AffineConstraint) 

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

365 

366 # Reset options. 

367 self.int.resetParams() 

368 

369 # verbosity 

370 picosVerbosity = self.verbosity() 

371 if picosVerbosity <= -1: 

372 scipVerbosity = 0 

373 elif picosVerbosity == 0: 

374 scipVerbosity = 2 

375 elif picosVerbosity == 1: 

376 scipVerbosity = 3 

377 elif picosVerbosity >= 2: 

378 scipVerbosity = 5 

379 self.int.setIntParam("display/verblevel", scipVerbosity) 

380 

381 # TODO: "numerics/lpfeastol" has been replaced in SCIP 7.0.0 (PySCIPOpt 

382 # 3.0.0) by "numerics/lpfeastolfactor", which is multiplied with 

383 # some unknown default feasibility (maybe "numerics/feastol"?). 

384 # Looking at it now it is not clear to me which of SCIP's 

385 # feasibility tolerances does what exactly and whether they are 

386 # absolute or relative measures. The workaround is to maintain 

387 # PICOS' old guess but ignore abs_prim_fsb_tol with the new SCIP. 

388 if int(scip.__version__.split(".")[0]) < 3: 

389 # abs_prim_fsb_tol 

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

391 self.int.setRealParam("numerics/lpfeastol", 

392 self.ext.options.abs_prim_fsb_tol) 

393 

394 # abs_dual_fsb_tol 

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

396 self.int.setRealParam("numerics/dualfeastol", 

397 self.ext.options.abs_dual_fsb_tol) 

398 

399 # rel_prim_fsb_tol, rel_dual_fsb_tol 

400 relFsbTols = [tol for tol in (self.ext.options.rel_prim_fsb_tol, 

401 self.ext.options.rel_dual_fsb_tol) if tol is not None] 

402 if relFsbTols: 

403 self.int.setRealParam("numerics/feastol", min(relFsbTols)) 

404 

405 # abs_ipm_opt_tol, abs_bnb_opt_tol 

406 absOptTols = [tol for tol in (self.ext.options.abs_ipm_opt_tol, 

407 self.ext.options.abs_bnb_opt_tol) if tol is not None] 

408 if absOptTols: 

409 self.int.setRealParam("limits/absgap", min(absOptTols)) 

410 

411 # rel_ipm_opt_tol, rel_bnb_opt_tol 

412 relOptTols = [tol for tol in (self.ext.options.rel_ipm_opt_tol, 

413 self.ext.options.rel_bnb_opt_tol) if tol is not None] 

414 if relOptTols: 

415 self.int.setRealParam("limits/gap", min(relOptTols)) 

416 

417 # timelimit 

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

419 self.int.setRealParam("limits/time", self.ext.options.timelimit) 

420 

421 # treememory 

422 if self.ext.options.treememory is not None: 

423 self.int.setRealParam("limits/memory", self.ext.options.treememory) 

424 

425 # max_fsb_nodes 

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

427 self.int.setRealParam( 

428 "limits/solutions", float(self.ext.options.max_fsb_nodes)) 

429 

430 # Handle SCIP-specific options. 

431 for key, value in self.ext.options.scip_params.items(): 

432 try: 

433 if isinstance(value, bool): 

434 self.int.setBoolParam(key, value) 

435 elif isinstance(value, str): 

436 if len(value) == 1: 

437 try: 

438 self.int.setCharParam(key, ord(value)) 

439 except LookupError: 

440 self.int.setStringParam(key, value) 

441 else: 

442 self.int.setStringParam(key, value) 

443 elif isinstance(value, float): 

444 self.int.setRealParam(key, value) 

445 elif isinstance(value, int): 

446 try: 

447 self.int.setIntParam(key, value) 

448 except LookupError: 

449 try: 

450 self.int.setLongintParam(key, value) 

451 except LookupError: 

452 self.int.setRealParam(key, float(value)) 

453 else: 

454 raise TypeError( 

455 "The value '{}' has an uncommon type.".format(value)) 

456 except KeyError as error: 

457 self._handle_bad_solver_specific_option_key(key, error) 

458 except ValueError as error: 

459 self._handle_bad_solver_specific_option_value(key, value, error) 

460 except (TypeError, LookupError) as error: 

461 picos_option = "{}_params".format(self.name) 

462 assert picos_option in self.ext.options, \ 

463 "The PICOS option '{}' does not exist.".format(picos_option) 

464 

465 raise OptionValueError( 

466 "Failed to guess type of {} option '{}' set via '{}'." 

467 .format(self.display_name, key, picos_option)) from error 

468 

469 # Handle unsupported options. 

470 self._handle_unsupported_options( 

471 "hotstart", "lp_node_method", "lp_root_method", "max_iterations") 

472 

473 # In the case of a pure LP, disable presolve to get duals. 

474 if self.ext.options.duals is not False and isLP: 

475 self._debug("Disabling SCIP's presolve, heuristics, and propagation" 

476 " features to allow for LP duals.") 

477 

478 self.int.setPresolve(scip.SCIP_PARAMSETTING.OFF) 

479 self.int.setHeuristics(scip.SCIP_PARAMSETTING.OFF) 

480 

481 # Note that this is a helper to set options, so they get reset at 

482 # the beginning of the function instead of in the else-scope below. 

483 self.int.disablePropagation() 

484 else: 

485 self.int.setPresolve(scip.SCIP_PARAMSETTING.DEFAULT) 

486 self.int.setHeuristics(scip.SCIP_PARAMSETTING.DEFAULT) 

487 

488 # Attempt to solve the problem. 

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

490 self.int.optimize() 

491 

492 # Retrieve primals. 

493 primals = {} 

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

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

496 value = [] 

497 for localIndex in range(picosVar.dim): 

498 scipVar = self._scipVar[picosVar.id_at(localIndex)] 

499 value.append(self.int.getVal(scipVar)) 

500 

501 primals[picosVar] = value 

502 

503 # Retrieve duals for LPs. 

504 duals = {} 

505 if self.ext.options.duals is not False and isLP: 

506 for picosConstraint in self.ext.constraints.values(): 

507 if isinstance(picosConstraint, DummyConstraint): 

508 duals[picosConstraint] = cvxopt.spmatrix( 

509 [], [], [], picosConstraint.size) 

510 continue 

511 

512 assert isinstance(picosConstraint, AffineConstraint) 

513 

514 # Retrieve dual value for constraint. 

515 scipDuals = [] 

516 for scipCon in self._scipCons[picosConstraint]: 

517 scipDuals.append(self.int.getDualsolLinear(scipCon)) 

518 picosDual = cvxopt.matrix(scipDuals, picosConstraint.size) 

519 

520 # Flip sign based on constraint relation. 

521 if not picosConstraint.is_increasing(): 

522 picosDual = -picosDual 

523 

524 # Flip sign based on objective sense. 

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

526 picosDual = -picosDual 

527 

528 duals[picosConstraint] = picosDual 

529 

530 # Retrieve objective value. 

531 value = self.int.getObjVal() 

532 

533 # Retrieve solution status. 

534 status = self.int.getStatus() 

535 if status == "optimal": 

536 primalStatus = SS_OPTIMAL 

537 dualStatus = SS_OPTIMAL 

538 problemStatus = PS_FEASIBLE 

539 elif status == "timelimit": 

540 primalStatus = SS_PREMATURE 

541 dualStatus = SS_PREMATURE 

542 problemStatus = PS_UNKNOWN 

543 elif status == "infeasible": 

544 primalStatus = SS_INFEASIBLE 

545 dualStatus = SS_UNKNOWN 

546 problemStatus = PS_INFEASIBLE 

547 elif status == "unbounded": 

548 primalStatus = SS_UNKNOWN 

549 dualStatus = SS_INFEASIBLE 

550 problemStatus = PS_UNBOUNDED 

551 else: 

552 primalStatus = SS_UNKNOWN 

553 dualStatus = SS_UNKNOWN 

554 problemStatus = PS_UNKNOWN 

555 

556 return self._make_solution( 

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

558 

559 

560# -------------------------------------- 

561__all__ = api_end(_API_START, globals())