Coverage for picos/solvers/solver_scip.py: 68.47%

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

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

20 

21import cvxopt 

22 

23from .. import settings 

24from ..apidoc import api_end, api_start 

25from ..constraints import (AffineConstraint, ConvexQuadraticConstraint, 

26 DummyConstraint, NonconvexQuadraticConstraint, 

27 RSOCConstraint, SOCConstraint) 

28from ..expressions import AffineExpression, BinaryVariable, IntegerVariable 

29from ..modeling.footprint import Specification 

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

31 PS_UNKNOWN, SS_INFEASIBLE, SS_OPTIMAL, 

32 SS_PREMATURE, SS_UNKNOWN) 

33from .solver import OptionValueError, Solver 

34 

35_API_START = api_start(globals()) 

36# ------------------------------- 

37 

38 

39class SCIPSolver(Solver): 

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

41 

42 SUPPORTED = Specification( 

43 objectives=[ 

44 AffineExpression], 

45 constraints=[ 

46 DummyConstraint, 

47 AffineConstraint, 

48 SOCConstraint, 

49 RSOCConstraint, 

50 ConvexQuadraticConstraint, 

51 NonconvexQuadraticConstraint]) 

52 

53 @classmethod 

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

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

56 result = Solver.supports(footprint, explain) 

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

58 return result 

59 

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

61 # produces dual values for LPs (which are sometimes incorrect). 

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

63 # with them at this point. 

64 if not settings.UNRELIABLE_STRATEGIES \ 

65 and footprint.continuous \ 

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

67 if explain: 

68 return (False, 

69 "Continuous problems, except for nonconvex quadratic ones " 

70 "(PICOS' choice).") 

71 else: 

72 return False 

73 

74 if footprint not in cls.SUPPORTED: 

75 if explain: 

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

77 else: 

78 return False 

79 

80 return (True, None) if explain else True 

81 

82 @classmethod 

83 def default_penalty(cls): 

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

85 return 1.0 # Stable free/open source solver. 

86 

87 @classmethod 

88 def test_availability(cls): 

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

90 cls.check_import("pyscipopt") 

91 

92 @classmethod 

93 def names(cls): 

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

95 return "scip", "SCIP", "SCIP Optimization Suite", "PySCIPOpt" 

96 

97 @classmethod 

98 def is_free(cls): 

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

100 return False 

101 

102 def __init__(self, problem): 

103 """Initialize a SCIP solver interface. 

104 

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

106 """ 

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

108 

109 self._scipVars = [] 

110 """A list of all SCIP variables added.""" 

111 

112 self._scipVarOffset = dict() 

113 """Maps PICOS variables to SCIP variable offsets.""" 

114 

115 self._scipCons = dict() 

116 """ 

117 Maps PICOS constraints to lists of SCIP constraints. 

118 

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

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

121 auxiliary constraint. 

122 """ 

123 

124 self._scipQuadObjAuxVar = None 

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

126 

127 self._scipQuadObjAuxCon = None 

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

129 

130 def reset_problem(self): 

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

132 self.int = None 

133 

134 self._scipVars.clear() 

135 self._scipVarOffset.clear() 

136 self._scipCons.clear() 

137 

138 self._scipQuadObjAuxVar = None 

139 self._scipQuadObjAuxCon = None 

140 

141 def _import_variable(self, picosVar): 

142 from math import ceil, floor 

143 

144 dim = picosVar.dim 

145 inf = self.int.infinity() 

146 

147 # Retrieve type code. 

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

149 scipVarType = "I" 

150 else: 

151 scipVarType = "C" 

152 

153 # Retrieve bounds. 

154 lowerBounds = [-inf]*dim 

155 upperBounds = [inf]*dim 

156 lower, upper = picosVar.bound_dicts 

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

158 lowerBounds[i] = b 

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

160 upperBounds[i] = b 

161 

162 # Refine bounds. 

163 if isinstance(picosVar, IntegerVariable): 

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

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

166 

167 # Import scalar variables. 

168 self._scipVarOffset[picosVar] = len(self._scipVars) 

169 for i in range(dim): 

170 self._scipVars.append(self.int.addVar( 

171 vtype=scipVarType, lb=lowerBounds[i], ub=upperBounds[i])) 

172 

173 def _import_variable_values(self): 

174 # TODO: Import variable values with SCIP. 

175 raise NotImplementedError 

176 

177 def _reset_variable_values(self): 

178 # TODO: Import variable values with SCIP. 

179 raise NotImplementedError 

180 

181 def _affinexp_pic2scip(self, picosExpression): 

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

183 

184 Multidimensional PICOS expressions are converted to a number of scalar 

185 SCIP expressions. 

186 

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

188 <pyscipopt.scip.Expr>`. 

189 """ 

190 from pyscipopt.scip import Expr 

191 

192 if picosExpression is None: 

193 yield Expr() 

194 return 

195 

196 rows = picosExpression.sparse_rows(self._scipVarOffset) 

197 for scipVars, coefficients, constant in rows: 

198 scipExpression = Expr() 

199 for scipVarIndex, coefficient in zip(scipVars, coefficients): 

200 scipExpression += coefficient * self._scipVars[scipVarIndex] 

201 scipExpression += constant 

202 yield scipExpression 

203 

204 def _scalar_affinexp_pic2scip(self, picosExpression): 

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

206 

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

208 """ 

209 assert len(picosExpression) == 1 

210 return next(self._affinexp_pic2scip(picosExpression)) 

211 

212 def _quadexp_pic2scip(self, picosExpression): 

213 # Convert affine part. 

214 scipExpression = self._scalar_affinexp_pic2scip(picosExpression.aff) 

215 

216 # Convert quadratic part. 

217 for (pVar1, pVar2), pCoefficients \ 

218 in picosExpression._sparse_quads.items(): 

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

220 localVar1Index = pCoefficients.I[sparseIndex] 

221 localVar2Index = pCoefficients.J[sparseIndex] 

222 localCoefficient = pCoefficients.V[sparseIndex] 

223 scipVar1 = self._scipVars[ 

224 self._scipVarOffset[pVar1] + localVar1Index] 

225 scipVar2 = self._scipVars[ 

226 self._scipVarOffset[pVar2] + localVar2Index] 

227 scipExpression += localCoefficient * scipVar1 * scipVar2 

228 

229 return scipExpression 

230 

231 def _import_linear_constraint(self, picosConstraint): 

232 assert isinstance(picosConstraint, AffineConstraint) 

233 

234 scipCons = [] 

235 picosLHS = picosConstraint.lmr 

236 for scipLHS in self._affinexp_pic2scip(picosLHS): 

237 if picosConstraint.is_increasing(): 

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

239 elif picosConstraint.is_decreasing(): 

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

241 elif picosConstraint.is_equality(): 

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

243 else: 

244 assert False, "Unexpected constraint relation." 

245 

246 return scipCons 

247 

248 def _import_quad_constraint(self, picosConstraint): 

249 assert isinstance(picosConstraint, 

250 (ConvexQuadraticConstraint, NonconvexQuadraticConstraint)) 

251 

252 scipLHS = self._quadexp_pic2scip(picosConstraint.le0) 

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

254 

255 return scipCons 

256 

257 def _import_socone_constraint(self, picosConstraint): 

258 scipCons = [] 

259 scipQuadLHS = self._quadexp_pic2scip( 

260 (picosConstraint.ne | picosConstraint.ne) 

261 - (picosConstraint.ub * picosConstraint.ub)) 

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

263 

264 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub) 

265 if scipAuxLHS.degree() > 0: 

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

267 

268 return scipCons 

269 

270 def _import_rscone_constraint(self, picosConstraint): 

271 scipCons = [] 

272 scipLHS = self._quadexp_pic2scip( 

273 (picosConstraint.ne | picosConstraint.ne) 

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

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

276 

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

278 # expressions is non-negative. 

279 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub1) 

280 if scipAuxLHS.degree() > 0: 

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

282 else: 

283 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub2) 

284 if scipAuxLHS.degree() > 0: 

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

286 

287 return scipCons 

288 

289 def _import_constraint(self, picosConstraint): 

290 # Import constraint based on type. 

291 if isinstance(picosConstraint, AffineConstraint): 

292 scipCons = self._import_linear_constraint(picosConstraint) 

293 elif isinstance(picosConstraint, 

294 (ConvexQuadraticConstraint, NonconvexQuadraticConstraint)): 

295 scipCons = self._import_quad_constraint(picosConstraint) 

296 elif isinstance(picosConstraint, SOCConstraint): 

297 scipCons = self._import_socone_constraint(picosConstraint) 

298 elif isinstance(picosConstraint, RSOCConstraint): 

299 scipCons = self._import_rscone_constraint(picosConstraint) 

300 else: 

301 assert isinstance(picosConstraint, DummyConstraint), \ 

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

303 picosConstraint.__class__.__name__) 

304 

305 # Map PICOS constraints to lists of SCIP constraints. 

306 self._scipCons[picosConstraint] = scipCons 

307 

308 def _import_objective(self): 

309 picosSense, picosObjective = self.ext.no 

310 

311 # Retrieve objective sense. 

312 if picosSense == "min": 

313 scipSense = "minimize" 

314 else: 

315 assert picosSense == "max" 

316 scipSense = "maximize" 

317 

318 # Import objective function. 

319 scipObjective = self._scalar_affinexp_pic2scip(picosObjective) 

320 

321 self.int.setObjective(scipObjective, scipSense) 

322 

323 def _import_problem(self): 

324 import pyscipopt as scip 

325 

326 # Create a problem instance. 

327 self.int = scip.Model() 

328 

329 # Import variables. 

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

331 self._import_variable(variable) 

332 

333 # Import constraints. 

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

335 self._import_constraint(constraint) 

336 

337 # Set objective. 

338 self._import_objective() 

339 

340 def _update_problem(self): 

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

342 raise NotImplementedError 

343 

344 def _solve(self): 

345 import pyscipopt as scip 

346 

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

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

349 and all(isinstance(constraint, AffineConstraint) 

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

351 

352 # Reset options. 

353 self.int.resetParams() 

354 

355 # verbosity 

356 picosVerbosity = self.verbosity() 

357 if picosVerbosity <= -1: 

358 scipVerbosity = 0 

359 elif picosVerbosity == 0: 

360 scipVerbosity = 2 

361 elif picosVerbosity == 1: 

362 scipVerbosity = 3 

363 elif picosVerbosity >= 2: 

364 scipVerbosity = 5 

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

366 

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

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

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

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

371 # feasibility tolerances does what exactly and whether they are 

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

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

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

375 # abs_prim_fsb_tol 

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

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

378 self.ext.options.abs_prim_fsb_tol) 

379 

380 # abs_dual_fsb_tol 

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

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

383 self.ext.options.abs_dual_fsb_tol) 

384 

385 # rel_prim_fsb_tol, rel_dual_fsb_tol 

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

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

388 if relFsbTols: 

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

390 

391 # abs_ipm_opt_tol, abs_bnb_opt_tol 

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

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

394 if absOptTols: 

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

396 

397 # rel_ipm_opt_tol, rel_bnb_opt_tol 

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

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

400 if relOptTols: 

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

402 

403 # timelimit 

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

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

406 

407 # treememory 

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

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

410 

411 # max_fsb_nodes 

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

413 self.int.setRealParam( 

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

415 

416 # Handle SCIP-specific options. 

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

418 try: 

419 if isinstance(value, bool): 

420 self.int.setBoolParam(key, value) 

421 elif isinstance(value, str): 

422 if len(value) == 1: 

423 try: 

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

425 except LookupError: 

426 self.int.setStringParam(key, value) 

427 else: 

428 self.int.setStringParam(key, value) 

429 elif isinstance(value, float): 

430 self.int.setRealParam(key, value) 

431 elif isinstance(value, int): 

432 try: 

433 self.int.setIntParam(key, value) 

434 except LookupError: 

435 try: 

436 self.int.setLongintParam(key, value) 

437 except LookupError: 

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

439 else: 

440 raise TypeError( 

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

442 except KeyError as error: 

443 self._handle_bad_solver_specific_option_key(key, error) 

444 except ValueError as error: 

445 self._handle_bad_solver_specific_option_value(key, value, error) 

446 except (TypeError, LookupError) as error: 

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

448 assert picos_option in self.ext.options, \ 

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

450 

451 raise OptionValueError( 

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

453 .format(self.short_name, key, picos_option)) from error 

454 

455 # Handle unsupported options. 

456 self._handle_unsupported_options( 

457 "hotstart", "lp_node_method", "lp_root_method", "max_iterations") 

458 

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

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

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

462 " features to allow for LP duals.") 

463 

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

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

466 

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

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

469 self.int.disablePropagation() 

470 else: 

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

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

473 

474 # Attempt to solve the problem. 

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

476 self.int.optimize() 

477 

478 # Retrieve primals. 

479 primals = {} 

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

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

482 scipVarOffset = self._scipVarOffset[picosVar] 

483 

484 value = [] 

485 for localIndex in range(picosVar.dim): 

486 scipVar = self._scipVars[scipVarOffset + localIndex] 

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

488 

489 primals[picosVar] = value 

490 

491 # Retrieve duals for LPs. 

492 duals = {} 

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

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

495 if isinstance(picosConstraint, DummyConstraint): 

496 duals[picosConstraint] = cvxopt.spmatrix( 

497 [], [], [], picosConstraint.size) 

498 continue 

499 

500 assert isinstance(picosConstraint, AffineConstraint) 

501 

502 # Retrieve dual value for constraint. 

503 scipDuals = [] 

504 for scipCon in self._scipCons[picosConstraint]: 

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

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

507 

508 # Flip sign based on constraint relation. 

509 if not picosConstraint.is_increasing(): 

510 picosDual = -picosDual 

511 

512 # Flip sign based on objective sense. 

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

514 picosDual = -picosDual 

515 

516 duals[picosConstraint] = picosDual 

517 

518 # Retrieve objective value. 

519 value = self.int.getObjVal() 

520 

521 # Retrieve solution status. 

522 status = self.int.getStatus() 

523 if status == "optimal": 

524 primalStatus = SS_OPTIMAL 

525 dualStatus = SS_OPTIMAL 

526 problemStatus = PS_FEASIBLE 

527 elif status == "timelimit": 

528 primalStatus = SS_PREMATURE 

529 dualStatus = SS_PREMATURE 

530 problemStatus = PS_UNKNOWN 

531 elif status == "infeasible": 

532 primalStatus = SS_INFEASIBLE 

533 dualStatus = SS_UNKNOWN 

534 problemStatus = PS_INFEASIBLE 

535 elif status == "unbounded": 

536 primalStatus = SS_UNKNOWN 

537 dualStatus = SS_INFEASIBLE 

538 problemStatus = PS_UNBOUNDED 

539 else: 

540 primalStatus = SS_UNKNOWN 

541 dualStatus = SS_UNKNOWN 

542 problemStatus = PS_UNKNOWN 

543 

544 return self._make_solution( 

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

546 

547 

548# -------------------------------------- 

549__all__ = api_end(_API_START, globals())