Coverage for picos/solvers/solver_gurobi.py: 77.55%
548 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +0000
« 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:`GurobiSolver`."""
21from collections import namedtuple
23import cvxopt
24import numpy
26from .. import settings
27from ..apidoc import api_end, api_start
28from ..constraints import (AffineConstraint, ConvexQuadraticConstraint,
29 DummyConstraint, NonconvexQuadraticConstraint,
30 RSOCConstraint, SOCConstraint)
31from ..expressions import (CONTINUOUS_VARTYPES, AffineExpression,
32 BinaryVariable, IntegerVariable,
33 QuadraticExpression)
34from ..expressions.data import cvx2csr, cvx2np
35from ..modeling.footprint import Specification
36from ..modeling.solution import (PS_FEASIBLE, PS_INF_OR_UNB, PS_INFEASIBLE,
37 PS_UNBOUNDED, PS_UNKNOWN, PS_UNSTABLE,
38 SS_EMPTY, SS_FEASIBLE, SS_INFEASIBLE,
39 SS_OPTIMAL, SS_PREMATURE, SS_UNKNOWN)
40from .solver import Solver
42_API_START = api_start(globals())
43# -------------------------------
46class GurobiSolver(Solver):
47 """Interface to the Gurobi solver via its official Python interface."""
49 # TODO: Don't support (conic) quadratic constraints when duals are
50 # requested because their precision is bad and can't be controlled?
51 SUPPORTED_8 = Specification(
52 objectives=[
53 AffineExpression,
54 QuadraticExpression],
55 constraints=[
56 DummyConstraint,
57 AffineConstraint,
58 SOCConstraint,
59 RSOCConstraint,
60 ConvexQuadraticConstraint])
62 SUPPORTED_9 = Specification(
63 objectives=[
64 AffineExpression,
65 QuadraticExpression],
66 constraints=[
67 DummyConstraint,
68 AffineConstraint,
69 SOCConstraint,
70 RSOCConstraint,
71 NonconvexQuadraticConstraint])
73 @classmethod
74 def _gurobi9(cls):
75 try:
76 import gurobipy as gurobi
77 except ImportError:
78 # This method should be used only after test_availability confirmed
79 # that gurobipy is available, however that method does not actually
80 # perform an import and so it could still fail here due to a bad
81 # installation. In this case an exception will be raised when Gurobi
82 # is actually selected and we just return False as a dummy here.
83 return False
84 else:
85 return hasattr(gurobi, "MVar")
87 @classmethod
88 def supports(cls, footprint, explain=False):
89 """Implement :meth:`~.solver.Solver.supports`."""
90 result = Solver.supports(footprint, explain)
91 if not result or (explain and not result[0]):
92 return result
94 supported = cls.SUPPORTED_9 if cls._gurobi9() else cls.SUPPORTED_8
96 if footprint not in supported:
97 if explain:
98 return False, supported.mismatch_reason(footprint)
99 else:
100 return False
102 return (True, None) if explain else True
104 @classmethod
105 def default_penalty(cls):
106 """Implement :meth:`~.solver.Solver.default_penalty`."""
107 return 0.0 # Commercial solver.
109 @classmethod
110 def test_availability(cls):
111 """Implement :meth:`~.solver.Solver.test_availability`."""
112 cls.check_import("gurobipy")
114 @classmethod
115 def names(cls):
116 """Implement :meth:`~.solver.Solver.names`."""
117 return "gurobi", "Gurobi", "Gurobi Optimizer", None
119 @classmethod
120 def is_free(cls):
121 """Implement :meth:`~.solver.Solver.is_free`."""
122 return False
124 GurobiMetaConstraint = namedtuple(
125 "GurobiMetaConstraint", ("auxCons", "auxVars"))
127 def __init__(self, problem):
128 """Initialize a Gurobi solver interface.
130 :param ~picos.Problem problem: The problem to be solved.
131 """
132 super(GurobiSolver, self).__init__(problem)
134 self._matint_decision = None
136 self._gurobiVar = dict()
137 """Maps PICOS variable indices to Gurobi variables. (Matrix interf.)"""
139 self._gurobiVars = []
140 """A list of all Gurobi variables added. (Legacy interface.)"""
142 self._gurobiVarOffset = dict()
143 """Maps PICOS variables to Gurobi variable offsets. (Legacy interf.)"""
145 self._gurobiLinCon = dict()
146 """Maps a PICOS linear constraint to (a) Gurobi linear constraint(s)."""
148 self._gurobiQuadCon = dict()
149 """Maps a PICOS quadr. constraint to (a) Gurobi quadr. constraint(s)."""
151 self._gurobiConicCon = dict()
152 """Maps a PICOS quadr. constraint to its Gurobi representation."""
154 def _make_matint_decision(self):
155 default = settings.PREFER_GUROBI_MATRIX_INTERFACE
156 choice = self.ext.options.gurobi_matint
158 if not choice and (not default or choice is not None):
159 return False
161 if not self._gurobi9():
162 if choice:
163 raise RuntimeError(
164 "Gurobi's matrix interface should be used by user's choice "
165 "but Gurobi < 9 appears to be installed. More precisely, "
166 "gurobipy.Mvar is not available.")
167 else:
168 return False
170 try:
171 self.check_import("scipy.sparse")
172 except ModuleNotFoundError as error:
173 if choice:
174 raise RuntimeError(
175 "Gurobi's matrix interface should be used by user's choice "
176 "but this requires SciPy, which was not found.") from error
177 else:
178 return False
180 assert choice or (choice is None and default)
181 return True
183 @property
184 def matint(self):
185 """Whether Gurobi's matrix interface is in use."""
186 if self._matint_decision is None:
187 self._matint_decision = self._make_matint_decision()
189 decision = self._matint_decision
190 choice = self.ext.options.gurobi_matint
192 if (choice and not decision) \
193 or (decision and not choice and choice is not None):
194 raise NotImplementedError(
195 "The user's choice with respect to using Gurobi's matrix "
196 "interface has changed between solution attempts. This is not "
197 "supported. To re-load the problem with the other interface, "
198 "you must manually reset your problem's solution strategy.")
200 return decision
202 def reset_problem(self):
203 """Implement :meth:`~.solver.Solver.reset_problem`."""
204 self.int = None
206 self._gurobiVar.clear()
207 self._gurobiVars.clear()
208 self._gurobiVarOffset.clear()
210 self._gurobiLinCon.clear()
211 self._gurobiQuadCon.clear()
212 self._gurobiConicCon.clear()
214 def _import_variable(self, picosVar):
215 import gurobipy as gurobi
217 dim = picosVar.dim
219 # Retrieve types.
220 if isinstance(picosVar, CONTINUOUS_VARTYPES):
221 gurobiVarType = gurobi.GRB.CONTINUOUS
222 elif isinstance(picosVar, IntegerVariable):
223 gurobiVarType = gurobi.GRB.INTEGER
224 elif isinstance(picosVar, BinaryVariable):
225 gurobiVarType = gurobi.GRB.BINARY
226 else:
227 assert False, "Unexpected variable type."
229 # Retrieve bounds.
230 lowerBounds = [-gurobi.GRB.INFINITY]*dim
231 upperBounds = [gurobi.GRB.INFINITY]*dim
232 lower, upper = picosVar.bound_dicts
233 for i, b in lower.items():
234 lowerBounds[i] = b
235 for i, b in upper.items():
236 upperBounds[i] = b
238 # Import the variable.
239 if self.matint:
240 gurobiVar = self.int.addMVar(dim, lb=lowerBounds, ub=upperBounds,
241 vtype=gurobiVarType, name=picosVar.name)
243 self._gurobiVar[picosVar] = gurobiVar
244 else:
245 gurobiVarsDict = self.int.addVars(
246 dim, lb=lowerBounds, ub=upperBounds, vtype=gurobiVarType)
247 gurobiVars = [gurobiVarsDict[i] for i in range(dim)]
249 self._gurobiVarOffset[picosVar] = len(self._gurobiVars)
250 self._gurobiVars.extend(gurobiVars)
252 def _remove_variable(self, picosVar):
253 if self.matint:
254 gurobiVar = self._gurobiVar.pop(picosVar)
256 self.int.remove(gurobiVar)
257 else:
258 offset = self._gurobiVarOffset[picosVar]
259 dim = picosVar.dim
261 gurobiVars = self._gurobiVars[offset:offset + dim]
263 self._gurobiVars = (
264 self._gurobiVars[:offset] + self._gurobiVars[offset + dim:])
266 for other in self._gurobiVarOffset:
267 if self._gurobiVarOffset[other] > offset:
268 self._gurobiVarOffset[other] -= dim
270 self.int.remove(gurobiVars)
272 def _import_variable_values(self):
273 for picosVar in self.ext.variables.values():
274 if picosVar.valued:
275 value = picosVar.internal_value
277 if self.matint:
278 gurobiVar = self._gurobiVar[picosVar]
279 gurobiVar.Start = value
280 else:
281 offset = self._gurobiVarOffset[picosVar]
282 dim = picosVar.dim
284 gurobiVars = self._gurobiVars[offset, offset + picosVar.dim]
286 for localIndex in range(dim):
287 gurobiVars[localIndex].Start = value[localIndex]
289 def _reset_variable_values(self):
290 import gurobipy as gurobi
292 if self.matint:
293 gurobiVars = self._gurobiVar.values()
294 else:
295 gurobiVars = self._gurobiVars
297 for gurobiVar in gurobiVars:
298 gurobiVar.Start = gurobi.GRB.UNDEFINED
300 def _affexp_pic2grb_matint(self, picosExpression):
301 assert self.matint
303 # NOTE: Constant Gurobi matrix expressions don't exist; return thus
304 # constant PICOS expressions as NumPy arrays.
305 gurobiExpression = numpy.ravel(cvx2np(
306 picosExpression._constant_coef))
308 for picosVar, coef in picosExpression._linear_coefs.items():
309 A = cvx2csr(coef)
310 x = self._gurobiVar[picosVar]
312 # NOTE: Using __(r)matmul__ as PICOS supports Python 3.4 and the
313 # @-operator was implemented in Python 3.5.
314 gurobiExpression += x.__rmatmul__(A)
316 return gurobiExpression
318 def _affexp_pic2grb_legacy(self, picosExpression):
319 import gurobipy as gurobi
321 assert not self.matint
323 for J, V, c in picosExpression.sparse_rows(self._gurobiVarOffset):
324 gurobiVars = [self._gurobiVars[j] for j in J]
325 gurobiExpression = gurobi.LinExpr(V, gurobiVars)
326 gurobiExpression.addConstant(c)
328 yield gurobiExpression
330 def _scalar_affexp_pic2grb(self, picosExpression):
331 assert len(picosExpression) == 1
333 if self.matint:
334 gurobiExpression = self._affexp_pic2grb_matint(picosExpression)
336 if picosExpression.constant:
337 assert isinstance(gurobiExpression, numpy.ndarray)
338 assert gurobiExpression.shape == (1,)
340 return gurobiExpression[0]
341 else:
342 return gurobiExpression
343 else:
344 return next(self._affexp_pic2grb_legacy(picosExpression))
346 def _quadexp_pic2grb(self, picosExpression):
347 import gurobipy as gurobi
349 assert isinstance(picosExpression, QuadraticExpression)
351 if self.matint:
352 # Import affine part.
353 gurobiExpression = self._affexp_pic2grb_matint(picosExpression.aff)
355 # Import quadratic part.
356 for picosVars, coef in picosExpression._sparse_quads.items():
357 Q = cvx2csr(coef)
358 x = self._gurobiVar[picosVars[0]]
359 y = self._gurobiVar[picosVars[1]]
361 gurobiExpression += x.__matmul__(Q).__matmul__(y)
362 else:
363 # Import affine part.
364 gurobiExpression = gurobi.QuadExpr(
365 self._scalar_affexp_pic2grb(picosExpression.aff))
367 # Import quadratic part.
368 V, I, J = [], [], []
369 for (x, y), Q in picosExpression._sparse_quads.items():
370 xOffset = self._gurobiVarOffset[x]
371 yOffset = self._gurobiVarOffset[y]
373 V.extend(Q.V)
374 I.extend(self._gurobiVars[i] for i in Q.I + xOffset)
375 J.extend(self._gurobiVars[j] for j in Q.J + yOffset)
377 gurobiExpression.addTerms(V, I, J)
379 return gurobiExpression
381 def _import_linear_constraint(self, picosCon):
382 import gurobipy as gurobi
384 assert isinstance(picosCon, AffineConstraint)
386 if self.matint:
387 gurobiLHS = self._affexp_pic2grb_matint(picosCon.lhs)
388 gurobiRHS = self._affexp_pic2grb_matint(picosCon.rhs)
390 # HACK: Fallback to the legacy interface for constant constraints.
391 # NOTE: This happens to work with remove_constraint since
392 # gurobipy.Model.remove accepts both lists and constraints.
393 if isinstance(gurobiLHS, numpy.ndarray) \
394 and isinstance(gurobiRHS, numpy.ndarray):
395 if picosCon.is_increasing():
396 gurobiSense = gurobi.GRB.LESS_EQUAL
397 elif picosCon.is_decreasing():
398 gurobiSense = gurobi.GRB.GREATER_EQUAL
399 elif picosCon.is_equality():
400 gurobiSense = gurobi.GRB.EQUAL
401 else:
402 assert False, "Unexpected constraint relation."
404 return [self.int.addLConstr(a, gurobiSense, b)
405 for a, b in zip(gurobiLHS, gurobiRHS)]
407 # Construct the constraint.
408 if picosCon.is_increasing():
409 gurobiCon = gurobiLHS <= gurobiRHS
410 elif picosCon.is_decreasing():
411 gurobiCon = gurobiLHS >= gurobiRHS
412 elif picosCon.is_equality():
413 gurobiCon = gurobiLHS == gurobiRHS
414 else:
415 assert False, "Unexpected constraint relation."
417 # Add the constraint.
418 gurobiCon = self.int.addConstr(gurobiCon)
420 return gurobiCon
421 else:
422 # Retrieve sense.
423 if picosCon.is_increasing():
424 gurobiSense = gurobi.GRB.LESS_EQUAL
425 elif picosCon.is_decreasing():
426 gurobiSense = gurobi.GRB.GREATER_EQUAL
427 elif picosCon.is_equality():
428 gurobiSense = gurobi.GRB.EQUAL
429 else:
430 assert False, "Unexpected constraint relation."
432 # Append scalar constraints.
433 gurobiCons = [self.int.addLConstr(gurobiLHS, gurobiSense, 0.0)
434 for gurobiLHS in self._affexp_pic2grb_legacy(picosCon.lmr)]
436 return gurobiCons
438 def _import_quad_constraint(self, picosCon):
439 import gurobipy as gurobi
441 # NOTE: NonconvexQuadraticConstraint includes ConvexQuadraticConstraint.
442 assert isinstance(picosCon, NonconvexQuadraticConstraint)
444 if self.matint:
445 gurobiLE0 = self._quadexp_pic2grb(picosCon.le0)
446 gurobiCon = self.int.addConstr(gurobiLE0 <= 0)
447 else:
448 gurobiLHS = self._quadexp_pic2grb(picosCon.le0)
449 gurobiRHS = -gurobiLHS.getLinExpr().getConstant()
451 if gurobiRHS:
452 gurobiLHS.getLinExpr().addConstant(gurobiRHS)
454 gurobiCon = self.int.addQConstr(
455 gurobiLHS, gurobi.GRB.LESS_EQUAL, gurobiRHS)
457 return gurobiCon
459 def _import_socone_constraint(self, picosCon):
460 import gurobipy as gurobi
462 assert isinstance(picosCon, SOCConstraint)
464 n = len(picosCon.ne)
466 # Load defining expressions.
467 gurobiRHS = self._scalar_affexp_pic2grb(picosCon.ub)
469 if self.matint:
470 # Load defining expressions.
471 gurobiLHS = self._affexp_pic2grb_matint(picosCon.ne)
473 # Add auxiliary variables for both sides.
474 gurobiLHSVar = self.int.addMVar(n, lb=-gurobi.GRB.INFINITY)
475 gurobiRHSVar = self.int.addMVar(1)
477 # Add constraints to identify auxiliary variables with expressions.
478 gurobiLHSCon = self.int.addConstr(gurobiLHSVar == gurobiLHS)
479 gurobiRHSCon = self.int.addConstr(gurobiRHSVar == gurobiRHS)
481 # Add a quadratic constraint over the auxiliary variables that
482 # represents the PICOS second order cone constraint itself.
483 gurobiQuadLHS = gurobiLHSVar.__matmul__(gurobiLHSVar)
484 gurobiQuadRHS = gurobiRHSVar.__matmul__(gurobiRHSVar)
485 gurobiQuadCon = self.int.addConstr(gurobiQuadLHS <= gurobiQuadRHS)
487 # Collect auxiliary objects.
488 auxCons = [gurobiLHSCon, gurobiRHSCon, gurobiQuadCon]
489 auxVars = [gurobiLHSVar, gurobiRHSVar]
490 else:
491 # Load defining expressions.
492 gurobiLHS = self._affexp_pic2grb_legacy(picosCon.ne)
494 # Add auxiliary variables: One for every dimension of the left hand
495 # side of the PICOS constraint and one for its right hand side.
496 gurobiLHSVarsDict = self.int.addVars(
497 n, lb=-gurobi.GRB.INFINITY, ub=gurobi.GRB.INFINITY)
498 gurobiLHSVars = gurobiLHSVarsDict.values()
499 gurobiRHSVar = self.int.addVar(lb=0.0, ub=gurobi.GRB.INFINITY)
501 # Add constraints that identify the left hand side Gurobi auxiliary
502 # variables with entries of the PICOS left hand side expression.
503 gurobiLHSDict = dict(enumerate(gurobiLHS))
504 gurobiLHSConsDict = self.int.addConstrs(
505 gurobiLHSVarsDict[d] == gurobiLHSDict[d] for d in range(n))
506 gurobiLHSCons = gurobiLHSConsDict.values()
508 # Add a constraint that identifies the right hand side Gurobi
509 # auxiliary variable with the PICOS right hand side expression.
510 gurobiRHSCon = self.int.addLConstr(
511 gurobiRHSVar, gurobi.GRB.EQUAL, gurobiRHS)
513 # Add a quadratic constraint over the auxiliary variables that
514 # represents the PICOS second order cone constraint itself.
515 quadExpr = gurobi.QuadExpr()
516 quadExpr.addTerms([1.0] * n, gurobiLHSVars, gurobiLHSVars)
517 gurobiQuadCon = self.int.addQConstr(
518 quadExpr, gurobi.GRB.LESS_EQUAL, gurobiRHSVar * gurobiRHSVar)
520 # Collect auxiliary objects.
521 auxCons = list(gurobiLHSCons) + [gurobiRHSCon, gurobiQuadCon]
522 auxVars = list(gurobiLHSVars) + [gurobiRHSVar]
524 return self.GurobiMetaConstraint(auxCons=auxCons, auxVars=auxVars)
526 def _import_rscone_constraint(self, picosCon):
527 import gurobipy as gurobi
529 assert isinstance(picosCon, RSOCConstraint)
531 n = len(picosCon.ne)
533 # Load defining expressions.
534 gurobiRHS = (
535 self._scalar_affexp_pic2grb(picosCon.ub1),
536 self._scalar_affexp_pic2grb(picosCon.ub2))
538 if self.matint:
539 # Load defining expressions.
540 gurobiLHS = self._affexp_pic2grb_matint(picosCon.ne)
542 # Add auxiliary variables for both sides.
543 gurobiLHSVar = self.int.addMVar(n, lb=-gurobi.GRB.INFINITY)
544 gurobiRHSVars = (self.int.addMVar(1), self.int.addMVar(1))
546 # Add constraints to identify auxiliary variables with expressions.
547 gurobiLHSCon = self.int.addConstr(gurobiLHSVar == gurobiLHS)
548 gurobiRHSConsDict = self.int.addConstrs(
549 gurobiRHSVars[i] == gurobiRHS[i] for i in range(2))
550 gurobiRHSCons = gurobiRHSConsDict.values()
552 # Add a quadratic constraint over the auxiliary variables that
553 # represents the PICOS rotated second order cone constraint itself.
554 gurobiQuadLHS = gurobiLHSVar.__matmul__(gurobiLHSVar)
555 gurobiQuadRHS = gurobiRHSVars[0].__matmul__(gurobiRHSVars[1])
556 gurobiQuadCon = self.int.addConstr(gurobiQuadLHS <= gurobiQuadRHS)
558 # Collect auxiliary objects.
559 auxCons = [gurobiLHSCon] + list(gurobiRHSCons) + [gurobiQuadCon]
560 auxVars = [gurobiLHSVar] + list(gurobiRHSVars)
561 else:
562 # Load defining expressions.
563 gurobiLHS = self._affexp_pic2grb_legacy(picosCon.ne)
565 # Add auxiliary variables: One for every dimension of the left hand
566 # side of the PICOS constraint and one for its right hand side.
567 gurobiLHSVarsDict = self.int.addVars(
568 n, lb=-gurobi.GRB.INFINITY, ub=gurobi.GRB.INFINITY)
569 gurobiLHSVars = gurobiLHSVarsDict.values()
570 gurobiRHSVars = self.int.addVars(
571 2, lb=0.0, ub=gurobi.GRB.INFINITY).values()
573 # Add constraints that identify the left hand side Gurobi auxiliary
574 # variables with entries of the PICOS left hand side expression.
575 gurobiLHSDict = dict(enumerate(gurobiLHS))
576 gurobiLHSConsDict = self.int.addConstrs(
577 gurobiLHSVarsDict[d] == gurobiLHSDict[d] for d in range(n))
578 gurobiLHSCons = gurobiLHSConsDict.values()
580 # Add two constraints that identify the right hand side Gurobi
581 # auxiliary variables with the PICOS right hand side expressions.
582 gurobiRHSConsDict = self.int.addConstrs(
583 gurobiRHSVars[i] == gurobiRHS[i] for i in (0, 1))
584 gurobiRHSCons = gurobiRHSConsDict.values()
586 # Add a quadratic constraint over the auxiliary variables that
587 # represents the PICOS second order cone constraint itself.
588 quadExpr = gurobi.QuadExpr()
589 quadExpr.addTerms([1.0] * n, gurobiLHSVars, gurobiLHSVars)
590 gurobiQuadCon = self.int.addQConstr(quadExpr, gurobi.GRB.LESS_EQUAL,
591 gurobiRHSVars[0] * gurobiRHSVars[1])
593 # Collect auxiliary objects.
594 auxCons = (
595 list(gurobiLHSCons) + list(gurobiRHSCons) + [gurobiQuadCon])
596 auxVars = list(gurobiLHSVars) + list(gurobiRHSVars)
598 return self.GurobiMetaConstraint(auxCons=auxCons, auxVars=auxVars)
600 def _import_constraint(self, picosCon):
601 if isinstance(picosCon, AffineConstraint):
602 self._gurobiLinCon[picosCon] = \
603 self._import_linear_constraint(picosCon)
604 elif isinstance(picosCon, NonconvexQuadraticConstraint):
605 self._gurobiQuadCon[picosCon] = \
606 self._import_quad_constraint(picosCon)
607 elif isinstance(picosCon, SOCConstraint):
608 self._gurobiConicCon[picosCon] = \
609 self._import_socone_constraint(picosCon)
610 elif isinstance(picosCon, RSOCConstraint):
611 self._gurobiConicCon[picosCon] = \
612 self._import_rscone_constraint(picosCon)
613 else:
614 assert isinstance(picosCon, DummyConstraint), \
615 "Unexpected constraint type: {}".format(
616 picosCon.__class__.__name__)
618 def _remove_constraint(self, picosCon):
619 if isinstance(picosCon, AffineConstraint):
620 self.int.remove(self._gurobiLinCon.pop(picosCon))
621 elif isinstance(picosCon, NonconvexQuadraticConstraint):
622 self.int.remove(self._gurobiQuadCon.pop(picosCon))
623 elif isinstance(picosCon, (SOCConstraint, RSOCConstraint)):
624 metaCon = self._gurobiConicCon.pop(picosCon)
626 self.int.remove(metaCon.auxCons)
627 self.int.remove(metaCon.auxVars)
628 else:
629 assert isinstance(picosCon, DummyConstraint), \
630 "Unexpected constraint type: {}".format(
631 picosCon.__class__.__name__)
633 def _import_objective(self):
634 import gurobipy as gurobi
636 picosSense, picosObjective = self.ext.no
638 # Retrieve objective sense.
639 if picosSense == "min":
640 gurobiSense = gurobi.GRB.MINIMIZE
641 else:
642 assert picosSense == "max"
643 gurobiSense = gurobi.GRB.MAXIMIZE
645 # Retrieve objective function.
646 if isinstance(picosObjective, AffineExpression):
647 gurobiObjective = self._scalar_affexp_pic2grb(picosObjective)
648 else:
649 assert isinstance(picosObjective, QuadraticExpression)
650 gurobiObjective = self._quadexp_pic2grb(picosObjective)
652 self.int.setObjective(gurobiObjective, gurobiSense)
654 def _import_problem(self):
655 import gurobipy as gurobi
657 # Create a problem instance.
658 if self._license_warnings:
659 self.int = gurobi.Model()
660 else:
661 with self._enforced_verbosity():
662 self.int = gurobi.Model()
664 # Import variables.
665 for variable in self.ext.variables.values():
666 self._import_variable(variable)
668 # Import constraints.
669 for constraint in self.ext.constraints.values():
670 self._import_constraint(constraint)
672 # Set objective.
673 self._import_objective()
675 def _update_problem(self):
676 for oldConstraint in self._removed_constraints():
677 self._remove_constraint(oldConstraint)
679 for oldVariable in self._removed_variables():
680 self._remove_variable(oldVariable)
682 for newVariable in self._new_variables():
683 self._import_variable(newVariable)
685 for newConstraint in self._new_constraints():
686 self._import_constraint(newConstraint)
688 if self._objective_has_changed():
689 self._import_objective()
691 def _solve(self):
692 import gurobipy as gurobi
694 # Reset options.
695 # NOTE: OutputFlag = 0 prevents resetParams from printing to console.
696 self.int.Params.OutputFlag = 0
697 self.int.resetParams()
699 # verbosity
700 self.int.Params.OutputFlag = 1 if self.verbosity() > 0 else 0
702 # abs_prim_fsb_tol
703 if self.ext.options.abs_prim_fsb_tol is not None:
704 self.int.Params.FeasibilityTol = self.ext.options.abs_prim_fsb_tol
706 # abs_dual_fsb_tol
707 if self.ext.options.abs_dual_fsb_tol is not None:
708 self.int.Params.OptimalityTol = self.ext.options.abs_dual_fsb_tol
710 # rel_ipm_opt_tol
711 if self.ext.options.rel_ipm_opt_tol is not None:
712 self.int.Params.BarConvTol = self.ext.options.rel_ipm_opt_tol
714 # HACK: Work around low precision (conic) quadratic duals.
715 self.int.Params.BarQCPConvTol = \
716 0.01 * self.ext.options.rel_ipm_opt_tol
718 # abs_bnb_opt_tol
719 if self.ext.options.abs_bnb_opt_tol is not None:
720 self.int.Params.MIPGapAbs = self.ext.options.abs_bnb_opt_tol
722 # rel_bnb_opt_tol
723 if self.ext.options.rel_bnb_opt_tol is not None:
724 self.int.Params.MIPGap = self.ext.options.rel_bnb_opt_tol
726 # integrality_tol
727 if self.ext.options.integrality_tol is not None:
728 self.int.Params.IntFeasTol = self.ext.options.integrality_tol
730 # markowitz_tol
731 if self.ext.options.markowitz_tol is not None:
732 self.int.Params.MarkowitzTol = self.ext.options.markowitz_tol
734 # max_iterations
735 if self.ext.options.max_iterations is not None:
736 self.int.Params.BarIterLimit = self.ext.options.max_iterations
737 self.int.Params.IterationLimit = self.ext.options.max_iterations
739 _lpm = {"interior": 2, "psimplex": 0, "dsimplex": 1}
741 # lp_node_method
742 if self.ext.options.lp_node_method is not None:
743 value = self.ext.options.lp_node_method
744 assert value in _lpm, "Unexpected lp_node_method value."
745 self.int.Params.SiftMethod = _lpm[value]
747 # lp_root_method
748 if self.ext.options.lp_root_method is not None:
749 value = self.ext.options.lp_root_method
750 assert value in _lpm, "Unexpected lp_root_method value."
751 self.int.Params.Method = _lpm[value]
753 # timelimit
754 if self.ext.options.timelimit is not None:
755 self.int.Params.TimeLimit = self.ext.options.timelimit
757 # max_fsb_nodes
758 if self.ext.options.max_fsb_nodes is not None:
759 self.int.Params.SolutionLimit = self.ext.options.max_fsb_nodes
761 # hotstart
762 if self.ext.options.hotstart:
763 self._import_variable_values()
764 else:
765 self._reset_variable_values()
767 # Handle Gurobi-specific options.
768 for key, value in self.ext.options.gurobi_params.items():
769 if not self.int.getParamInfo(key):
770 self._handle_bad_solver_specific_option_key(key)
772 try:
773 self.int.setParam(key, value)
774 except TypeError as error:
775 self._handle_bad_solver_specific_option_value(key, value, error)
777 # Handle unsupported options.
778 self._handle_unsupported_option("treememory")
780 # Extend functionality for continuous problems.
781 if self.ext.is_continuous():
782 # Compute duals also for QPs and QC(Q)Ps.
783 if self.ext.options.duals is not False:
784 self.int.setParam(gurobi.GRB.Param.QCPDual, 1)
786 # Allow nonconvex quadratic objectives.
787 # TODO: Allow querying self.ext.objective directly.
788 # TODO: Check if this should/must be set also for Gurobi >= 9.
789 if self.ext.footprint.nonconvex_quadratic_objective:
790 self.int.setParam(gurobi.GRB.Param.NonConvex, 2)
792 # Attempt to solve the problem.
793 with self._header(), self._stopwatch():
794 try:
795 self.int.optimize()
796 except gurobi.GurobiError as error:
797 if error.errno == gurobi.GRB.Error.Q_NOT_PSD:
798 self._handle_continuous_nonconvex_error(error)
799 else:
800 raise
802 # Retrieve primals.
803 primals = {}
804 if self.ext.options.primals is not False:
805 for picosVar in self.ext.variables.values():
806 try:
807 if self.matint:
808 value = cvxopt.matrix(self._gurobiVar[picosVar].X)
809 else:
810 o = self._gurobiVarOffset[picosVar]
811 d = picosVar.dim
813 value = [v.X for v in self._gurobiVars[o:o + d]]
814 except (AttributeError, gurobi.GurobiError):
815 # NOTE: AttributeError is raised for gurobipy.Var,
816 # gurobi.GurobiError for gurobipy.MVar.
817 primals[picosVar] = None
818 else:
819 primals[picosVar] = value
821 # Retrieve duals.
822 duals = {}
823 if self.ext.options.duals is not False and self.ext.is_continuous():
824 for picosCon in self.ext.constraints.values():
825 if isinstance(picosCon, DummyConstraint):
826 duals[picosCon] = cvxopt.spmatrix([], [], [], picosCon.size)
827 continue
829 # HACK: Work around gurobiCon.getAttr(gurobi.GRB.Attr.Pi)
830 # printing a newline to console when it raises an
831 # AttributeError and OutputFlag is enabled. This is a
832 # WONTFIX on Gurobi's end (PICOS #264, Gurobi #14248).
833 # TODO: Check if this also happens for urobiCon.Pi, which is now
834 # used for both interfaces.
835 oldOutput = self.int.Params.OutputFlag
836 self.int.Params.OutputFlag = 0
838 try:
839 if isinstance(picosCon, AffineConstraint):
840 gurobiCon = self._gurobiLinCon[picosCon]
842 # HACK: Seee _import_linear_constraint.
843 if not self.matint or isinstance(gurobiCon, list):
844 gurobiDual = [c.Pi for c in gurobiCon]
845 else:
846 gurobiDual = gurobiCon.Pi
848 picosDual = cvxopt.matrix(gurobiDual, picosCon.size)
850 if not picosCon.is_increasing():
851 picosDual = -picosDual
852 elif isinstance(picosCon, SOCConstraint):
853 gurobiMetaCon = self._gurobiConicCon[picosCon]
855 if self.matint:
856 ne, ub, _ = gurobiMetaCon.auxCons
857 dual = numpy.hstack([ub.Pi, ne.Pi])
858 picosDual = cvxopt.matrix(dual)
859 else:
860 n = len(picosCon.ne)
861 assert len(gurobiMetaCon.auxCons) == n + 2
863 ne = gurobiMetaCon.auxCons[:n]
864 ub = gurobiMetaCon.auxCons[n]
866 z, lbd = [c.Pi for c in ne], ub.Pi
867 picosDual = cvxopt.matrix([lbd] + z)
868 elif isinstance(picosCon, RSOCConstraint):
869 gurobiMetaCon = self._gurobiConicCon[picosCon]
871 if self.matint:
872 ne, ub1, ub2, _ = gurobiMetaCon.auxCons
873 dual = numpy.hstack([ub1.Pi, ub2.Pi, ne.Pi])
874 picosDual = cvxopt.matrix(dual)
875 else:
876 n = len(picosCon.ne)
877 assert len(gurobiMetaCon.auxCons) == n + 3
879 ne = gurobiMetaCon.auxCons[:n]
880 ub1 = gurobiMetaCon.auxCons[n]
881 ub2 = gurobiMetaCon.auxCons[n + 1]
883 z, a, b = [c.Pi for c in ne], ub1.Pi, ub2.Pi
884 picosDual = cvxopt.matrix([a] + [b] + z)
885 elif isinstance(picosCon, NonconvexQuadraticConstraint):
886 picosDual = None
887 else:
888 assert isinstance(picosCon, DummyConstraint), \
889 "Unexpected constraint type: {}".format(
890 picosCon.__class__.__name__)
892 # Flip sign based on objective sense.
893 if picosDual and self.ext.no.direction == "min":
894 picosDual = -picosDual
895 except (AttributeError, gurobi.GurobiError):
896 # NOTE: AttributeError is raised for gurobipy.Constr,
897 # gurobi.GurobiError for gurobipy.MConstr.
898 duals[picosCon] = None
899 else:
900 duals[picosCon] = picosDual
902 # HACK: See above. Also: Silence Gurobi while enabling output.
903 if oldOutput != 0:
904 with self._enforced_verbosity(noStdOutAt=float("inf")):
905 self.int.Params.OutputFlag = oldOutput
907 # Retrieve objective value.
908 try:
909 value = self.int.ObjVal
910 except AttributeError:
911 value = None
913 # Retrieve solution status.
914 statusCode = self.int.Status
915 if statusCode == gurobi.GRB.Status.LOADED:
916 raise RuntimeError("Gurobi claims to have just loaded the problem "
917 "while PICOS expects the solution search to have terminated.")
918 elif statusCode == gurobi.GRB.Status.OPTIMAL:
919 primalStatus = SS_OPTIMAL
920 dualStatus = SS_OPTIMAL
921 problemStatus = PS_FEASIBLE
922 elif statusCode == gurobi.GRB.Status.INFEASIBLE:
923 primalStatus = SS_INFEASIBLE
924 dualStatus = SS_UNKNOWN
925 problemStatus = PS_INFEASIBLE
926 elif statusCode == gurobi.GRB.Status.INF_OR_UNBD:
927 primalStatus = SS_UNKNOWN
928 dualStatus = SS_UNKNOWN
929 problemStatus = PS_INF_OR_UNB
930 elif statusCode == gurobi.GRB.Status.UNBOUNDED:
931 primalStatus = SS_UNKNOWN
932 dualStatus = SS_INFEASIBLE
933 problemStatus = PS_UNBOUNDED
934 elif statusCode == gurobi.GRB.Status.CUTOFF:
935 # "Optimal objective for model was proven to be worse than the value
936 # specified in the Cutoff parameter. No solution information is
937 # available."
938 primalStatus = SS_PREMATURE
939 dualStatus = SS_PREMATURE
940 problemStatus = PS_UNKNOWN
941 elif statusCode == gurobi.GRB.Status.ITERATION_LIMIT:
942 primalStatus = SS_PREMATURE
943 dualStatus = SS_PREMATURE
944 problemStatus = PS_UNKNOWN
945 elif statusCode == gurobi.GRB.Status.NODE_LIMIT:
946 primalStatus = SS_PREMATURE
947 dualStatus = SS_EMPTY # Applies only to mixed integer problems.
948 problemStatus = PS_UNKNOWN
949 elif statusCode == gurobi.GRB.Status.TIME_LIMIT:
950 primalStatus = SS_PREMATURE
951 dualStatus = SS_PREMATURE
952 problemStatus = PS_UNKNOWN
953 elif statusCode == gurobi.GRB.Status.SOLUTION_LIMIT:
954 primalStatus = SS_PREMATURE
955 dualStatus = SS_PREMATURE
956 problemStatus = PS_UNKNOWN
957 elif statusCode == gurobi.GRB.Status.INTERRUPTED:
958 primalStatus = SS_PREMATURE
959 dualStatus = SS_PREMATURE
960 problemStatus = PS_UNKNOWN
961 elif statusCode == gurobi.GRB.Status.NUMERIC:
962 primalStatus = SS_UNKNOWN
963 dualStatus = SS_UNKNOWN
964 problemStatus = PS_UNSTABLE
965 elif statusCode == gurobi.GRB.Status.SUBOPTIMAL:
966 # "Unable to satisfy optimality tolerances; a sub-optimal solution
967 # is available."
968 primalStatus = SS_FEASIBLE
969 dualStatus = SS_FEASIBLE
970 problemStatus = PS_FEASIBLE
971 elif statusCode == gurobi.GRB.Status.INPROGRESS:
972 raise RuntimeError("Gurobi claims solution search to be 'in "
973 "progress' while PICOS expects it to have terminated.")
974 elif statusCode == gurobi.GRB.Status.USER_OBJ_LIMIT:
975 # "User specified an objective limit (a bound on either the best
976 # objective or the best bound), and that limit has been reached."
977 primalStatus = SS_FEASIBLE
978 dualStatus = SS_EMPTY # Applies only to mixed integer problems.
979 problemStatus = PS_FEASIBLE
980 else:
981 primalStatus = SS_UNKNOWN
982 dualStatus = SS_UNKNOWN
983 problemStatus = PS_UNKNOWN
985 return self._make_solution(
986 value, primals, duals, primalStatus, dualStatus, problemStatus)
989# --------------------------------------
990__all__ = api_end(_API_START, globals())