Coverage for picos/solvers/solver_cvxopt.py: 83.14%
338 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2012-2017 Guillaume Sagnol
3# Copyright (C) 2017-2022 Maximilian Stahlberg
4#
5# This file is part of PICOS.
6#
7# PICOS is free software: you can redistribute it and/or modify it under the
8# terms of the GNU General Public License as published by the Free Software
9# Foundation, either version 3 of the License, or (at your option) any later
10# version.
11#
12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
14# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License along with
17# this program. If not, see <http://www.gnu.org/licenses/>.
18# ------------------------------------------------------------------------------
20"""Implementation of :class:`CVXOPTSolver`."""
22import cvxopt
23import numpy
25from ..apidoc import api_end, api_start
26from ..constraints import (AffineConstraint, DummyConstraint, LMIConstraint,
27 LogSumExpConstraint, RSOCConstraint, SOCConstraint)
28from ..expressions import CONTINUOUS_VARTYPES, AffineExpression, LogSumExp
29from ..modeling.footprint import Specification
30from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED,
31 PS_UNKNOWN, SS_INFEASIBLE, SS_OPTIMAL,
32 SS_UNKNOWN, Solution)
33from .solver import Solver
35_API_START = api_start(globals())
36# -------------------------------
39class CVXOPTSolver(Solver):
40 """Interface to the CVXOPT solver.
42 Also used as an interface to the
43 :class:`SMCP solver <picos.solvers.solver_smcp.SMCPSolver>`.
44 """
46 SUPPORTED = Specification(
47 objectives=[
48 AffineExpression,
49 LogSumExp],
50 variables=CONTINUOUS_VARTYPES,
51 constraints=[
52 DummyConstraint,
53 AffineConstraint,
54 SOCConstraint,
55 RSOCConstraint,
56 LMIConstraint,
57 LogSumExpConstraint])
59 @classmethod
60 def supports(cls, footprint, explain=False):
61 """Implement :meth:`~.solver.Solver.supports`."""
62 result = Solver.supports(footprint, explain)
63 if not result or (explain and not result[0]):
64 return result
66 if footprint not in cls.SUPPORTED:
67 if explain:
68 return False, cls.SUPPORTED.mismatch_reason(footprint)
69 else:
70 return False
72 return (True, None) if explain else True
74 @classmethod
75 def default_penalty(cls):
76 """Implement :meth:`~.solver.Solver.default_penalty`."""
77 return 1.0 # Stable free/open source solver.
79 @classmethod
80 def test_availability(cls):
81 """Implement :meth:`~.solver.Solver.test_availability`."""
82 # CVXOPT is a dependency of PICOS, so it is always available.
83 pass
85 @classmethod
86 def names(cls):
87 """Implement :meth:`~.solver.Solver.names`."""
88 return "cvxopt", "CVXOPT", "Python Convex Optimization Solver", None
90 @classmethod
91 def is_free(cls):
92 """Implement :meth:`~.solver.Solver.is_free`."""
93 return True
95 @property
96 def is_smcp(self):
97 """Whether to implement SMCP instead of CVXOPT."""
98 return False
100 def __init__(self, problem):
101 """Initialize a CVXOPT solver interface.
103 :param ~picos.Problem problem: The problem to be solved.
104 """
105 super(CVXOPTSolver, self).__init__(problem)
107 self._numVars = 0
108 """Total number of scalar variables passed to CVXOPT."""
110 self._cvxoptVarOffset = {}
111 """Maps a PICOS variable to its offset in the constraint matrix."""
113 # HACK: Setting this to false would result in variable bounds to be
114 # ignored, instead of added to the linear inequalities matrix.
115 # This is used by a Problem method that prints the problem to a
116 # file using a CVXOPTSolver instance.
117 self.import_variable_bounds = True
119 # SMCP currently offers no option reset, so keep a backup.
120 if self.is_smcp:
121 import smcp
122 self._smcp_default_options = smcp.solvers.options.copy()
124 def reset_problem(self):
125 """Implement :meth:`~.solver.Solver.reset_problem`."""
126 self.int = None
128 self._numVars = 0
129 self._cvxoptVarOffset.clear()
131 def _is_geometric_program(self):
132 return isinstance(self.ext.no.function, LogSumExp) or \
133 any(isinstance(constraint, LogSumExpConstraint)
134 for constraint in self.ext.constraints.values())
136 def _affine_expression_to_G_and_h(self, expression):
137 assert isinstance(expression, AffineExpression)
139 return expression.sparse_matrix_form(
140 varOffsetMap=self._cvxoptVarOffset, dense_b=True)
142 _Gh = _affine_expression_to_G_and_h
144 def _import_variables_without_bounds(self):
145 offset = 0
146 for variable in self.ext.variables.values():
147 dim = variable.dim
149 # Register the variable.
150 self._cvxoptVarOffset[variable] = offset
151 offset += dim
153 assert offset == self._numVars
155 def _import_variable_bounds(self):
156 for variable in self.ext.variables.values():
157 if self.import_variable_bounds:
158 bounds = variable.bound_constraint
159 if bounds:
160 self._import_affine_constraint(bounds)
162 def _import_objective(self):
163 direction, objective = self.ext.no
165 if self._is_geometric_program():
166 if isinstance(objective, LogSumExp):
167 (F, g) = self._Gh(objective.x)
168 else:
169 assert isinstance(objective, AffineExpression)
170 (F, g) = self._Gh(objective)
172 # NOTE: Needs to happen before LogSumExpConstraint are added.
173 self.int["K"] = [F.size[0]]
175 if direction == "max":
176 F, g = -F, -g
178 self.int["F"] = F
179 self.int["g"] = g
180 else:
181 (c, _) = self._Gh(objective)
183 if direction == "max":
184 c = -c
186 # Must be a dense column-vector.
187 c = cvxopt.matrix(c).T
189 self.int["c"] = c
191 def _import_affine_constraint(self, constraint):
192 assert isinstance(constraint, AffineConstraint)
194 (G_smaller, h_smaller) = self._Gh(constraint.smaller)
195 (G_greater, h_greater) = self._Gh(constraint.greater)
197 G = G_smaller - G_greater
198 h = h_greater - h_smaller
200 if constraint.is_equality():
201 self.int["A"] = cvxopt.sparse([self.int["A"], G])
202 self.int["b"] = cvxopt.matrix([self.int["b"], h])
203 else:
204 self.int["Gl"] = cvxopt.sparse([self.int["Gl"], G])
205 self.int["hl"] = cvxopt.matrix([self.int["hl"], h])
207 def _import_soc_constraint(self, constraint):
208 assert isinstance(constraint, SOCConstraint)
210 (A, b) = self._Gh(constraint.ne)
211 (c, d) = self._Gh(constraint.ub)
213 self.int["Gq"].append(cvxopt.sparse([-c, -A]))
214 self.int["hq"].append(cvxopt.matrix([d, b]))
216 def _import_rsoc_constraint(self, constraint):
217 assert isinstance(constraint, RSOCConstraint)
219 (A, b) = self._Gh(constraint.ne)
220 (c1, d1) = self._Gh(constraint.ub1)
221 (c2, d2) = self._Gh(constraint.ub2)
223 self.int["Gq"].append(cvxopt.sparse([-c1 - c2, -2 * A, c2 - c1]))
224 self.int["hq"].append(cvxopt.matrix([d1 + d2, 2 * b, d1 - d2]))
226 def _import_lmi_constraint(self, constraint):
227 assert isinstance(constraint, LMIConstraint)
229 (G_smaller, h_smaller) = self._Gh(constraint.smaller)
230 (G_greater, h_greater) = self._Gh(constraint.greater)
232 self.int["Gs"].append(G_smaller - G_greater)
233 self.int["hs"].append(h_greater - h_smaller)
235 def _import_lse_constraint(self, constraint):
236 assert isinstance(constraint, LogSumExpConstraint)
238 (F, g) = self._Gh(constraint.le0.x)
240 self.int["F"] = cvxopt.sparse([self.int["F"], F])
241 self.int["g"] = cvxopt.matrix([self.int["g"], g])
242 self.int["K"].append(F.size[0])
244 def _import_constraint(self, constraint):
245 if isinstance(constraint, AffineConstraint):
246 self._import_affine_constraint(constraint)
247 elif isinstance(constraint, SOCConstraint):
248 self._import_soc_constraint(constraint)
249 elif isinstance(constraint, RSOCConstraint):
250 self._import_rsoc_constraint(constraint)
251 elif isinstance(constraint, LMIConstraint):
252 self._import_lmi_constraint(constraint)
253 elif isinstance(constraint, LogSumExpConstraint):
254 self._import_lse_constraint(constraint)
255 else:
256 assert isinstance(constraint, DummyConstraint), \
257 "Unexpected constraint type: {}".format(
258 constraint.__class__.__name__)
260 def _import_problem(self):
261 self._numVars = sum(var.dim for var in self.ext.variables.values())
263 # CVXOPT's internal problem representation is stateless; a number of
264 # matrices are supplied to the appropriate solver function each time a
265 # search is started. These matrices are thus stored in self.int.
266 self.int = {
267 # Objective function coefficients.
268 "c": None,
270 # Linear equality left hand side.
271 "A": cvxopt.spmatrix([], [], [], (0, self._numVars), tc="d"),
273 # Linear equality right hand side.
274 "b": cvxopt.matrix([], (0, 1), tc="d"),
276 # Linear inequality left hand side.
277 "Gl": cvxopt.spmatrix([], [], [], (0, self._numVars), tc="d"),
279 # Linear inequality right hand side.
280 "hl": cvxopt.matrix([], (0, 1), tc="d"),
282 # Second order cone inequalities left hand sides.
283 "Gq": [],
285 # Second order cone inequalities right hand sides.
286 "hq": [],
288 # Semidefinite cone inequalities left hand sides.
289 "Gs": [],
291 # Semidefinite cone inequalities right hand sides.
292 "hs": [],
294 # Geometric program data.
295 "F": None,
296 "g": None,
297 "K": None
298 }
300 # Import variables without their bounds.
301 self._import_variables_without_bounds()
303 # Set objective.
304 # NOTE: This needs to happen before constraints are added as
305 # self.int["K"][0] refers to an LSE objective while
306 # self.int["K"][i] with i > 0 refers to LSE constraints.
307 self._import_objective()
309 # Import constraints.
310 for constraint in self.ext.constraints.values():
311 self._import_constraint(constraint)
313 # Import variable bounds as additional affine constraints.
314 # NOTE: This needs to happen after constraints are added due to how the
315 # dual values for constraints are extracted.
316 self._import_variable_bounds()
318 def _update_problem(self):
319 raise NotImplementedError
321 def _solve(self):
322 if self.is_smcp:
323 import smcp
325 p = self.int
326 isGP = self._is_geometric_program()
328 # Clear all options set previously. This is necessary because CVXOPT
329 # options are global, and might be changed even by another problem.
330 cvxopt.solvers.options.clear()
332 # verbosity
333 cvxopt.solvers.options["show_progress"] = (self.verbosity() >= 1)
335 # rel_prim_fsb_tol, rel_dual_fsb_tol
336 feasibilityTols = [tol for tol in (self.ext.options.rel_prim_fsb_tol,
337 self.ext.options.rel_dual_fsb_tol) if tol is not None]
338 if feasibilityTols:
339 cvxopt.solvers.options["feastol"] = min(feasibilityTols)
341 # abs_ipm_opt_tol
342 if self.ext.options.abs_ipm_opt_tol is not None:
343 cvxopt.solvers.options["abstol"] = self.ext.options.abs_ipm_opt_tol
345 # rel_ipm_opt_tol
346 if self.ext.options.rel_ipm_opt_tol is not None:
347 cvxopt.solvers.options["reltol"] = self.ext.options.rel_ipm_opt_tol
349 # max_iterations
350 if self.ext.options.max_iterations is not None:
351 cvxopt.solvers.options["maxiters"] = self.ext.options.max_iterations
352 else:
353 cvxopt.solvers.options["maxiters"] = int(1e6)
355 # cvxopt_kktsolver
356 if self.ext.options.cvxopt_kktsolver is not None:
357 userKKT = self.ext.options.cvxopt_kktsolver
358 else:
359 userKKT = None
361 # cvxopt_kktreg
362 if self.ext.options.cvxopt_kktreg is not None:
363 cvxopt.solvers.options["kktreg"] = self.ext.options.cvxopt_kktreg
365 # Handle unsupported options.
366 self._handle_unsupported_options(
367 "lp_root_method", "lp_node_method", "timelimit", "treememory",
368 "max_fsb_nodes", "hotstart")
370 # TODO: Add CVXOPT-sepcific options. Candidates are:
371 # - refinement
373 # Set options for SMCP.
374 if self.is_smcp:
375 # Restore default options.
376 smcp.solvers.options = self._smcp_default_options.copy()
378 # Copy options also used by CVXOPT.
379 smcp.solvers.options.update(cvxopt.solvers.options)
381 # Further handle "verbose" option.
382 smcp.solvers.options["debug"] = (self.verbosity() >= 2)
384 # TODO: Add SMCP-sepcific options.
386 if self._debug():
387 from pprint import pformat
388 self._debug("Setting options:\n{}\n".format(pformat(
389 smcp.solvers.options if self.is_smcp
390 else cvxopt.solvers.options)))
392 # Print a header.
393 if self.is_smcp:
394 subsolverText = None
395 else:
396 if isGP:
397 subsolverText = "internal GP solver"
398 else:
399 subsolverText = "internal CONELP solver"
401 # Further prepare the problem for the CVXOPT/SMCP CONELP solvers.
402 # TODO: This should be done during import.
403 if not isGP:
404 # Retrieve the structure of the cone, which is a cartesian product
405 # of the non-negative orthant of dimension l, a number of second
406 # order cones with dimensions in q and a number of positive
407 # semidefinite cones with dimensions in s.
408 dims = {
409 "l": p["Gl"].size[0],
410 "q": [Gqi.size[0] for Gqi in p["Gq"]],
411 "s": [int(numpy.sqrt(Gsi.size[0])) for Gsi in p["Gs"]]
412 }
414 # Construct G and h to contain all conic inequalities, starting with
415 # those with respect to the non-negative orthant.
416 G = p["Gl"]
417 h = p["hl"]
419 # SMCP's ConeLP solver does not handle (linear) equalities, so cast
420 # them as inequalities.
421 if self.is_smcp:
422 smcp_eps = min(feasibilityTols)
423 if p["A"].size[0] > 0:
424 G = cvxopt.sparse([G, p["A"]])
425 G = cvxopt.sparse([G, -p["A"]])
426 h = cvxopt.matrix([h, p["b"]+smcp_eps])
427 h = cvxopt.matrix([h, smcp_eps-p["b"]])
428 dims["l"] += (2 * p["A"].size[0])
430 # Remove the lines in G and h corresponding to 0==0 or 0<=0
431 JP = list(set(G.I))
432 IP = range(len(JP))
433 VP = [1] * len(JP)
435 if len(JP) != dims["l"]:
436 # is there a constraint of the form 0<=a, (a<0) ?
437 if any([b < -smcp_eps for (i, b) in enumerate(h)
438 if i not in JP]):
439 raise Exception(
440 'infeasible constraint of the form '
441 '0 <= a, with a<0')
443 # left-multiply with PPP-matrix to remove 0-constraints
444 PPP = cvxopt.spmatrix(VP, IP, JP, (len(IP), G.size[0]))
445 dims["l"] = len(JP)
446 G = PPP * G
447 h = PPP * h
449 # Add second-order cone inequalities.
450 for i in range(len(dims["q"])):
451 G = cvxopt.sparse([G, p["Gq"][i]])
452 h = cvxopt.matrix([h, p["hq"][i]])
454 # Add semidefinite cone inequalities.
455 for i in range(len(dims["s"])):
456 G = cvxopt.sparse([G, p["Gs"][i]])
457 h = cvxopt.matrix([h, p["hs"][i]])
459 # Remove zero lines from linear equality constraint matrix, as
460 # CVXOPT expects this matrix to have full row rank.
461 JP = list(set(p["A"].I))
462 IP = range(len(JP))
463 VP = [1] * len(JP)
465 # Skip solution on an infeasible constraint.
466 if any([b for (i, b) in enumerate(p["b"]) if i not in JP]):
467 return Solution(
468 primals=None, duals=None, problem=self.ext, solver="PICOS",
469 primalStatus=SS_INFEASIBLE, dualStatus=SS_UNKNOWN,
470 problemStatus=PS_INFEASIBLE, vectorizedPrimals=True)
472 P = cvxopt.spmatrix(VP, IP, JP, (len(IP), p["A"].size[0]))
473 A = P * p["A"]
474 b = P * p["b"]
476 # Attempt to solve the problem.
477 with self._header(subsolverText), self._stopwatch():
478 if self.is_smcp:
479 if self._debug():
480 self._debug("Calling smcp.solvers.conelp(c, G, h, dims) "
481 "with\nc:\n{}\nG:\n{}\nh:\n{}\ndims:\n{}\n"
482 .format(p["c"], G, h, dims))
483 try:
484 result = smcp.solvers.conelp(p["c"], G, h, dims)
485 except TypeError:
486 # HACK: Work around "'NoneType' object is not subscriptable"
487 # exception with infeasible/unbounded problems.
488 result = None
489 else:
490 kwargs = {}
492 if userKKT:
493 kwargs["kktsolver"] = userKKT
494 elif isGP:
495 # Use the more reliable LDL solver right away.
496 kwargs["kktsolver"] = "ldl"
497 else:
498 # Try the fast but unreliable CHOL solver first.
499 kwargs["kktsolver"] = "chol"
501 if isGP:
502 result = cvxopt.solvers.gp(p["K"], p["F"], p["g"], p["Gl"],
503 p["hl"], p["A"], p["b"], **kwargs)
504 else:
505 try:
506 result = cvxopt.solvers.conelp(
507 p["c"], G, h, dims, A, b, **kwargs)
509 if not userKKT and result["status"] == "unknown":
510 raise ValueError("The first solution attempt with "
511 "CHOL as a KKT solver returnd a solution with "
512 "unknown status. This exception triggers "
513 "another solution attempt using LDL.")
514 except ValueError:
515 # NOTE: Apart from the one created by PICOS above, there
516 # are at least two more ValueError produced by
517 # CVXOPT when an unreliable KKT solver is used.
518 if userKKT:
519 raise # Always respect the user's choice.
520 else:
521 # Re-solve using the LDL solver.
522 # TODO: Consider pre-solving on PICOS end to prevent
523 # the "Rank(A) < p or Rank([G; A]) < n" error.
524 kwargs["kktsolver"] = "ldl"
525 result = cvxopt.solvers.conelp(
526 p["c"], G, h, dims, A, b, **kwargs)
528 # Retrieve primals.
529 primals = {}
530 if self.ext.options.primals is not False and result is not None \
531 and result["x"] is not None:
532 for variable in self.ext.variables.values():
533 offset = self._cvxoptVarOffset[variable]
534 value = list(result["x"][offset:offset + variable.dim])
535 primals[variable] = value
537 # Retrieve duals.
538 duals = {}
539 if self.ext.options.duals is not False and result is not None:
540 (indy, indzl, indzq, indznl, indzs) = (0, 0, 0, 0, 0)
542 if isGP:
543 zkey = "zl"
544 zqkey = "zq"
545 zskey = "zs"
546 else:
547 zkey = "z"
548 zqkey = "z"
549 zskey = "z"
550 indzq = dims["l"]
551 indzs = dims["l"] + sum(dims["q"])
553 if self.is_smcp:
554 # Equality constraints were cast as two inequalities.
555 ieq = p["Gl"].size[0]
556 neq = (dims["l"] - ieq) // 2
557 soleq = result["z"][ieq:ieq + neq]
558 soleq -= result["z"][ieq + neq:ieq + 2 * neq]
559 else:
560 soleq = result["y"]
562 for constraint in self.ext.constraints.values():
563 if isinstance(constraint, DummyConstraint):
564 duals[constraint] = cvxopt.spmatrix(
565 [], [], [], constraint.size)
566 continue
568 dual = None
569 consSz = len(constraint)
571 if isinstance(constraint, AffineConstraint):
572 if constraint.is_equality():
573 if soleq is not None:
574 dual = -(P.T * soleq)[indy:indy + consSz]
575 indy += consSz
576 else:
577 if result[zkey] is not None:
578 dual = result[zkey][indzl:indzl + consSz]
579 indzl += consSz
580 elif isinstance(constraint, SOCConstraint) \
581 or isinstance(constraint, RSOCConstraint):
582 if result[zqkey] is not None:
583 if isGP:
584 dual = result[zqkey][indzq]
585 dual[1:] = -dual[1:]
586 indzq += 1
587 else:
588 dual = result[zqkey][indzq:indzq + consSz]
589 if isinstance(constraint, RSOCConstraint):
590 # RScone were cast as a SOcone on import, so
591 # transform the dual to a proper RScone dual.
592 alpha = dual[0] + dual[-1]
593 beta = dual[0] - dual[-1]
594 z = 2.0 * dual[1:-1]
595 dual = cvxopt.matrix([alpha, beta, z])
596 indzq += consSz
597 elif isinstance(constraint, LMIConstraint):
598 if result[zskey] is not None:
599 matSz = constraint.size[0]
600 if isGP:
601 dual = cvxopt.matrix(
602 result[zskey][indzs], (matSz, matSz))
603 indzs += 1
604 else:
605 dual = cvxopt.matrix(
606 result[zskey][indzs:indzs + consSz],
607 (matSz, matSz))
608 indzs += consSz
609 elif isinstance(constraint, LogSumExpConstraint):
610 # TODO: Retrieve LSE duals.
611 indznl += 1
613 duals[constraint] = dual
615 # Retrieve objective value.
616 if result is None:
617 value = None
618 elif isGP:
619 value = None
620 else:
621 p = result['primal objective']
622 d = result['dual objective']
624 if p is not None and d is not None:
625 value = 0.5 * (p + d)
626 elif p is not None:
627 value = p
628 elif d is not None:
629 value = d
630 else:
631 value = None
633 if value is not None and self.ext.no.direction == "max":
634 value = -value
636 if self.is_smcp:
637 value = -value
639 # Retrieve solution status.
640 status = result["status"] if result else "unknown"
641 if status == "optimal":
642 primalStatus = SS_OPTIMAL
643 dualStatus = SS_OPTIMAL
644 problemStatus = PS_FEASIBLE
645 elif status == "primal infeasible":
646 primalStatus = SS_INFEASIBLE
647 dualStatus = SS_UNKNOWN
648 problemStatus = PS_INFEASIBLE
649 elif status == "dual infeasible":
650 primalStatus = SS_UNKNOWN
651 dualStatus = SS_INFEASIBLE
652 problemStatus = PS_UNBOUNDED
653 else:
654 primalStatus = SS_UNKNOWN
655 dualStatus = SS_UNKNOWN
656 problemStatus = PS_UNKNOWN
658 return self._make_solution(value, primals, duals, primalStatus,
659 dualStatus, problemStatus, {"cvxopt_sol": result})
662# --------------------------------------
663__all__ = api_end(_API_START, globals())