Coverage for picos/solvers/solver_scip.py: 68.47%
295 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2018-2022 Maximilian Stahlberg
3#
4# This file is part of PICOS.
5#
6# PICOS is free software: you can redistribute it and/or modify it under the
7# terms of the GNU General Public License as published by the Free Software
8# Foundation, either version 3 of the License, or (at your option) any later
9# version.
10#
11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
17# ------------------------------------------------------------------------------
19"""Implementation of :class:`SCIPSolver`."""
21import cvxopt
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
35_API_START = api_start(globals())
36# -------------------------------
39class SCIPSolver(Solver):
40 """Interface to the SCIP solver via the PySCIPOpt Python interface."""
42 SUPPORTED = Specification(
43 objectives=[
44 AffineExpression],
45 constraints=[
46 DummyConstraint,
47 AffineConstraint,
48 SOCConstraint,
49 RSOCConstraint,
50 ConvexQuadraticConstraint,
51 NonconvexQuadraticConstraint])
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
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
74 if footprint not in cls.SUPPORTED:
75 if explain:
76 return False, cls.SUPPORTED.mismatch_reason(footprint)
77 else:
78 return False
80 return (True, None) if explain else True
82 @classmethod
83 def default_penalty(cls):
84 """Implement :meth:`~.solver.Solver.default_penalty`."""
85 return 1.0 # Stable free/open source solver.
87 @classmethod
88 def test_availability(cls):
89 """Implement :meth:`~.solver.Solver.test_availability`."""
90 cls.check_import("pyscipopt")
92 @classmethod
93 def names(cls):
94 """Implement :meth:`~.solver.Solver.names`."""
95 return "scip", "SCIP", "SCIP Optimization Suite", "PySCIPOpt"
97 @classmethod
98 def is_free(cls):
99 """Implement :meth:`~.solver.Solver.is_free`."""
100 return False
102 def __init__(self, problem):
103 """Initialize a SCIP solver interface.
105 :param ~picos.Problem problem: The problem to be solved.
106 """
107 super(SCIPSolver, self).__init__(problem)
109 self._scipVars = []
110 """A list of all SCIP variables added."""
112 self._scipVarOffset = dict()
113 """Maps PICOS variables to SCIP variable offsets."""
115 self._scipCons = dict()
116 """
117 Maps PICOS constraints to lists of SCIP constraints.
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 """
124 self._scipQuadObjAuxVar = None
125 """A SCIP auxiliary variable to support a PICOS quadratic objective."""
127 self._scipQuadObjAuxCon = None
128 """A SCIP auxiliary constraint to support a PICOS quadr. objective."""
130 def reset_problem(self):
131 """Implement :meth:`~.solver.Solver.reset_problem`."""
132 self.int = None
134 self._scipVars.clear()
135 self._scipVarOffset.clear()
136 self._scipCons.clear()
138 self._scipQuadObjAuxVar = None
139 self._scipQuadObjAuxCon = None
141 def _import_variable(self, picosVar):
142 from math import ceil, floor
144 dim = picosVar.dim
145 inf = self.int.infinity()
147 # Retrieve type code.
148 if isinstance(picosVar, (BinaryVariable, IntegerVariable)):
149 scipVarType = "I"
150 else:
151 scipVarType = "C"
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
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]
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]))
173 def _import_variable_values(self):
174 # TODO: Import variable values with SCIP.
175 raise NotImplementedError
177 def _reset_variable_values(self):
178 # TODO: Import variable values with SCIP.
179 raise NotImplementedError
181 def _affinexp_pic2scip(self, picosExpression):
182 """Transform an affine expression from PICOS to SCIP.
184 Multidimensional PICOS expressions are converted to a number of scalar
185 SCIP expressions.
187 :returns: A :class:`list` of :class:`SCIP expressions
188 <pyscipopt.scip.Expr>`.
189 """
190 from pyscipopt.scip import Expr
192 if picosExpression is None:
193 yield Expr()
194 return
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
204 def _scalar_affinexp_pic2scip(self, picosExpression):
205 """Transform a PICOS scalar affine expression into a SCIP expression.
207 :returns: A :class:`SCIP expression <pyscipopt.scip.Expr>`.
208 """
209 assert len(picosExpression) == 1
210 return next(self._affinexp_pic2scip(picosExpression))
212 def _quadexp_pic2scip(self, picosExpression):
213 # Convert affine part.
214 scipExpression = self._scalar_affinexp_pic2scip(picosExpression.aff)
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
229 return scipExpression
231 def _import_linear_constraint(self, picosConstraint):
232 assert isinstance(picosConstraint, AffineConstraint)
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."
246 return scipCons
248 def _import_quad_constraint(self, picosConstraint):
249 assert isinstance(picosConstraint,
250 (ConvexQuadraticConstraint, NonconvexQuadraticConstraint))
252 scipLHS = self._quadexp_pic2scip(picosConstraint.le0)
253 scipCons = [self.int.addCons(scipLHS <= 0.0)]
255 return scipCons
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))
264 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub)
265 if scipAuxLHS.degree() > 0:
266 scipCons.append(self.int.addCons(scipAuxLHS >= 0.0))
268 return scipCons
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))
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))
287 return scipCons
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__)
305 # Map PICOS constraints to lists of SCIP constraints.
306 self._scipCons[picosConstraint] = scipCons
308 def _import_objective(self):
309 picosSense, picosObjective = self.ext.no
311 # Retrieve objective sense.
312 if picosSense == "min":
313 scipSense = "minimize"
314 else:
315 assert picosSense == "max"
316 scipSense = "maximize"
318 # Import objective function.
319 scipObjective = self._scalar_affinexp_pic2scip(picosObjective)
321 self.int.setObjective(scipObjective, scipSense)
323 def _import_problem(self):
324 import pyscipopt as scip
326 # Create a problem instance.
327 self.int = scip.Model()
329 # Import variables.
330 for variable in self.ext.variables.values():
331 self._import_variable(variable)
333 # Import constraints.
334 for constraint in self.ext.constraints.values():
335 self._import_constraint(constraint)
337 # Set objective.
338 self._import_objective()
340 def _update_problem(self):
341 # TODO: Support all problem updates supported by SCIP.
342 raise NotImplementedError
344 def _solve(self):
345 import pyscipopt as scip
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())
352 # Reset options.
353 self.int.resetParams()
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)
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)
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)
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))
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))
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))
403 # timelimit
404 if self.ext.options.timelimit is not None:
405 self.int.setRealParam("limits/time", self.ext.options.timelimit)
407 # treememory
408 if self.ext.options.treememory is not None:
409 self.int.setRealParam("limits/memory", self.ext.options.treememory)
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))
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)
451 raise OptionValueError(
452 "Failed to guess type of {} option '{}' set via '{}'."
453 .format(self.short_name, key, picos_option)) from error
455 # Handle unsupported options.
456 self._handle_unsupported_options(
457 "hotstart", "lp_node_method", "lp_root_method", "max_iterations")
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.")
464 self.int.setPresolve(scip.SCIP_PARAMSETTING.OFF)
465 self.int.setHeuristics(scip.SCIP_PARAMSETTING.OFF)
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)
474 # Attempt to solve the problem.
475 with self._header(), self._stopwatch():
476 self.int.optimize()
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]
484 value = []
485 for localIndex in range(picosVar.dim):
486 scipVar = self._scipVars[scipVarOffset + localIndex]
487 value.append(self.int.getVal(scipVar))
489 primals[picosVar] = value
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
500 assert isinstance(picosConstraint, AffineConstraint)
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)
508 # Flip sign based on constraint relation.
509 if not picosConstraint.is_increasing():
510 picosDual = -picosDual
512 # Flip sign based on objective sense.
513 if picosDual and self.ext.no.direction == "min":
514 picosDual = -picosDual
516 duals[picosConstraint] = picosDual
518 # Retrieve objective value.
519 value = self.int.getObjVal()
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
544 return self._make_solution(
545 value, primals, duals, primalStatus, dualStatus, problemStatus)
548# --------------------------------------
549__all__ = api_end(_API_START, globals())