Coverage for picos/solvers/solver_ecos.py: 83.78%
296 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:`ECOSSolver`."""
21import cvxopt
22import numpy
24from ..apidoc import api_end, api_start
25from ..constraints import (AffineConstraint, DummyConstraint,
26 ExpConeConstraint, 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, PS_UNSTABLE, SS_FAILURE,
31 SS_INFEASIBLE, SS_OPTIMAL, SS_PREMATURE,
32 SS_UNKNOWN)
33from .solver import Solver
35_API_START = api_start(globals())
36# -------------------------------
39class ECOSSolver(Solver):
40 """Interface to the ECOS solver via its official Python interface."""
42 SUPPORTED = Specification(
43 objectives=[
44 AffineExpression],
45 constraints=[
46 DummyConstraint,
47 AffineConstraint,
48 SOCConstraint,
49 RSOCConstraint,
50 ExpConeConstraint])
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
59 if footprint not in cls.SUPPORTED:
60 if explain:
61 return False, cls.SUPPORTED.mismatch_reason(footprint)
62 else:
63 return False
65 return (True, None) if explain else True
67 @classmethod
68 def default_penalty(cls):
69 """Implement :meth:`~.solver.Solver.default_penalty`."""
70 return 1.0 # Stable free/open source solver.
72 @classmethod
73 def test_availability(cls):
74 """Implement :meth:`~.solver.Solver.test_availability`."""
75 cls.check_import("ecos")
77 @classmethod
78 def names(cls):
79 """Implement :meth:`~.solver.Solver.names`."""
80 return "ecos", "ECOS", "Embedded Conic Solver", "ecos-python"
82 @classmethod
83 def is_free(cls):
84 """Implement :meth:`~.solver.Solver.is_free`."""
85 return True
87 def __init__(self, problem):
88 """Initialize an ECOS solver interface.
90 :param ~picos.Problem problem: The problem to be solved.
91 """
92 super(ECOSSolver, self).__init__(problem)
94 self._numVars = 0
95 """Total number of scalar variables passed to ECOS."""
97 self._ecosVarOffset = {}
98 """Maps a PICOS variable to its offset in the ECOS constraint matrix."""
100 self._ecosConIndices = dict()
101 """Maps a PICOS constraints to a range of ECOS constraint indices.
103 Note that equality constraints use a different index space."""
105 self._ecosModuleCache = None
106 """A cached reference to the ECOS module."""
108 def reset_problem(self):
109 """Implement :meth:`~.solver.Solver.reset_problem`."""
110 self.int = None
112 self._numVars = 0
113 self._ecosVarOffset.clear()
114 self._ecosConIndices.clear()
116 @property
117 def ecos(self):
118 """Return the ECOS core module ecos.py.
120 The module is obtained by ``import ecos`` up to ECOS 2.0.6 and by
121 ``import ecos.ecos`` starting with ECOS 2.0.7.
122 """
123 if self._ecosModuleCache is None:
124 import ecos
125 if hasattr(ecos, "ecos"):
126 self._ecosModuleCache = ecos.ecos
127 else:
128 self._ecosModuleCache = ecos
129 return self._ecosModuleCache
131 @staticmethod
132 def stack(*args):
133 """Stack vectors or matrices, the latter vertically."""
134 import scipy.sparse
136 if isinstance(args[0], scipy.sparse.spmatrix):
137 for i in range(1, len(args)):
138 assert isinstance(args[i], scipy.sparse.spmatrix)
140 return scipy.sparse.vstack(args, format="csc")
141 else:
142 for i in range(len(args)):
143 assert isinstance(args[i], numpy.ndarray)
145 return numpy.concatenate(args)
147 def _affine_expression_to_G_and_h(self, expression):
148 assert isinstance(expression, AffineExpression)
150 return expression.scipy_sparse_matrix_form(
151 varOffsetMap=self._ecosVarOffset, dense_b=True)
153 _Gh = _affine_expression_to_G_and_h
155 def _import_variables(self):
156 offset = 0
157 for variable in self.ext.variables.values():
158 dim = variable.dim
160 # Mark integer variables.
161 if isinstance(variable, (BinaryVariable, IntegerVariable)):
162 ecosIndices = range(offset, offset + dim)
164 if isinstance(variable, BinaryVariable):
165 self.int["bool_vars_idx"].extend(ecosIndices)
166 else:
167 self.int["int_vars_idx"].extend(ecosIndices)
169 # Register the variable.
170 self._ecosVarOffset[variable] = offset
171 offset += dim
173 # Import bounds.
174 for variable in self.ext.variables.values():
175 bounds = variable.bound_constraint
176 if bounds:
177 self._import_affine_constraint(bounds)
179 assert offset == self._numVars
181 def _append_equality(self, A, b):
182 offset = len(self.int["b"])
183 length = len(b)
185 self.int["A"] = self.stack(self.int["A"], A)
186 self.int["b"] = self.stack(self.int["b"], b)
188 return range(offset, offset + length)
190 def _append_inequality(self, G, h, typecode):
191 assert typecode in self.int["dims"]
193 # Make sure constraints are appened in the proper order.
194 if typecode == "l":
195 assert len(self.int["dims"]["q"]) == 0 \
196 and self.int["dims"]["e"] == 0
197 elif typecode == "q":
198 assert self.int["dims"]["e"] == 0
200 offset = len(self.int["h"])
201 length = len(h)
203 self.int["G"] = self.stack(self.int["G"], G)
204 self.int["h"] = self.stack(self.int["h"], h)
206 if typecode == "q":
207 self.int["dims"][typecode].append(length)
208 elif typecode == "e":
209 self.int["dims"][typecode] += 1
210 else:
211 self.int["dims"][typecode] += length
213 return range(offset, offset + length)
215 def _import_affine_constraint(self, constraint):
216 assert isinstance(constraint, AffineConstraint)
218 G, minus_h = self._Gh(constraint.le0)
219 h = -minus_h
221 if constraint.is_equality():
222 return self._append_equality(G, h)
223 else:
224 return self._append_inequality(G, h, "l")
226 def _import_socone_constraint(self, constraint):
227 assert isinstance(constraint, SOCConstraint)
229 (A, b) = self._Gh(constraint.ne)
230 (c, d) = self._Gh(constraint.ub)
232 return self._append_inequality(
233 self.stack(-c, -A), self.stack(d, b), "q")
235 def _import_rscone_constraint(self, constraint):
236 assert isinstance(constraint, RSOCConstraint)
238 (A, b) = self._Gh(constraint.ne)
239 (c1, d1) = self._Gh(constraint.ub1)
240 (c2, d2) = self._Gh(constraint.ub2)
242 return self._append_inequality(
243 self.stack(-(c1 + c2), -2.0*A, c2 - c1),
244 self.stack(d1 + d2, 2.0*b, d1 - d2), "q")
246 def _import_expcone_constraint(self, constraint):
247 assert isinstance(constraint, ExpConeConstraint)
249 (Gx, hx) = self._Gh(constraint.x)
250 (Gy, hy) = self._Gh(constraint.y)
251 (Gz, hz) = self._Gh(constraint.z)
253 # ECOS' exponential cone is cl{(x,y,z) | exp(x/z) <= y/z, z > 0},
254 # PICOS' is cl{(x,y,z) | x >= y*exp(z/y), y > 0}. Note that given y > 0
255 # it is x >= y*exp(z/y) if and only if exp(z/y) <= x/y. Therefor we can
256 # transform from our coordinates to theirs with the mapping
257 # (x, y, z) ↦ (z, x, y). Further, G and h with G = (Gx, Gy, Gz) and
258 # h = (hx, hy, hz) are such that G*X + h = (x, y, z) where X is the
259 # row-vectorization of all variables. ECOS however expects G and h such
260 # that h - G*X is constrained to be in the exponential cone.
261 return self._append_inequality(
262 -self.stack(Gz, Gx, Gy), self.stack(hz, hx, hy), "e")
264 def _import_objective(self):
265 direction, objective = self.ext.no
267 # ECOS only supports minimization; flip the sign for maximization.
268 sign = 1.0 if direction == "min" else -1.0
270 # Import coefficients.
271 for variable, coefficient in objective._linear_coefs.items():
272 for localIndex in range(variable.dim):
273 index = self._ecosVarOffset[variable] + localIndex
274 self.int["c"][index] = sign * coefficient[localIndex]
276 def _import_problem(self):
277 import scipy.sparse
279 self._numVars = sum(var.dim for var in self.ext.variables.values())
281 # ECOS' internal problem representation is stateless; a number of
282 # vectors and matrices is supplied each time a search is started.
283 # These vectors and matrices are thus stored in self.int.
284 self.int = {
285 # Objective function coefficients.
286 "c": numpy.zeros(self._numVars),
288 # Linear equality left hand side.
289 "A": scipy.sparse.csc_matrix((0, self._numVars)),
291 # Linear equality right hand side.
292 "b": numpy.zeros(0),
294 # Conic inequality left hand side.
295 "G": scipy.sparse.csc_matrix((0, self._numVars)),
297 # Conic inequality right hand side.
298 "h": numpy.zeros(0),
300 # Cone definition: Linear, second order, exponential dimensions.
301 "dims": {"l": 0, "q": [], "e": 0},
303 # Boolean variable indices.
304 "bool_vars_idx": [],
306 # Integer variable indices.
307 "int_vars_idx": []
308 }
310 # NOTE: Constraints need to be append ordered by type.
311 # TODO: Add fast constraints-by-type iterators to Problem.
313 # Import variables with their bounds as affine constraints.
314 self._import_variables()
316 # Import affine constraints.
317 for constraint in self.ext.constraints.values():
318 if isinstance(constraint, AffineConstraint):
319 self._ecosConIndices[constraint] = \
320 self._import_affine_constraint(constraint)
322 # Import second order cone constraints.
323 for constraint in self.ext.constraints.values():
324 if isinstance(constraint, SOCConstraint):
325 self._ecosConIndices[constraint] = \
326 self._import_socone_constraint(constraint)
328 # Import rotated second order cone constraints.
329 for constraint in self.ext.constraints.values():
330 if isinstance(constraint, RSOCConstraint):
331 self._ecosConIndices[constraint] = \
332 self._import_rscone_constraint(constraint)
334 # Import exponential cone constraints.
335 for constraint in self.ext.constraints.values():
336 if isinstance(constraint, ExpConeConstraint):
337 self._ecosConIndices[constraint] = \
338 self._import_expcone_constraint(constraint)
340 # Make sure that no unsupported constraints are present.
341 for constraint in self.ext.constraints.values():
342 assert isinstance(constraint, (DummyConstraint, AffineConstraint,
343 SOCConstraint, RSOCConstraint, ExpConeConstraint))
345 # Set objective.
346 self._import_objective()
348 def _update_problem(self):
349 raise NotImplementedError
351 def _solve(self):
352 ecosOptions = {}
354 # verbosity
355 beVerbose = self.verbosity() >= 1
356 ecosOptions["verbose"] = beVerbose
357 ecosOptions["mi_verbose"] = beVerbose
359 # rel_prim_fsb_tol, rel_dual_fsb_tol
360 feasibilityTols = [tol for tol in (self.ext.options.rel_prim_fsb_tol,
361 self.ext.options.rel_dual_fsb_tol) if tol is not None]
362 if feasibilityTols:
363 ecosOptions["feastol"] = min(feasibilityTols)
365 # abs_ipm_opt_tol
366 if self.ext.options.abs_ipm_opt_tol is not None:
367 ecosOptions["abstol"] = self.ext.options.abs_ipm_opt_tol
369 # rel_ipm_opt_tol
370 if self.ext.options.rel_ipm_opt_tol is not None:
371 ecosOptions["reltol"] = self.ext.options.rel_ipm_opt_tol
373 # abs_bnb_opt_tol
374 if self.ext.options.abs_bnb_opt_tol is not None:
375 ecosOptions["mi_abs_eps"] = self.ext.options.abs_bnb_opt_tol
377 # rel_bnb_opt_tol
378 if self.ext.options.rel_bnb_opt_tol is not None:
379 ecosOptions["mi_rel_eps"] = self.ext.options.rel_bnb_opt_tol
381 # integrality_tol
382 # HACK: ECOS_BB uses ECOS' "tolerance for feasibility condition for
383 # inaccurate solution" as the integrality tolerance.
384 if self.ext.options.integrality_tol is not None:
385 ecosOptions["feastol_inacc"] = self.ext.options.integrality_tol
387 # max_iterations
388 if self.ext.options.max_iterations is not None:
389 ecosOptions["max_iters"] = self.ext.options.max_iterations
390 ecosOptions["mi_max_iters"] = self.ext.options.max_iterations
392 # Handle unsupported options.
393 self._handle_unsupported_options(
394 "lp_root_method", "lp_node_method", "timelimit", "treememory",
395 "max_fsb_nodes", "hotstart")
397 # Assemble the solver input.
398 arguments = {}
399 arguments.update(self.int)
400 arguments.update(ecosOptions)
402 # Debug print the solver input.
403 if self._debug():
404 from pprint import pformat
405 self._debug("Passing to ECOS:\n" + pformat(arguments))
407 # Attempt to solve the problem.
408 with self._header(), self._stopwatch():
409 solution = self.ecos.solve(**arguments)
411 # Debug print the solver output.
412 if self._debug():
413 from pprint import pformat
414 self._debug("Recevied from ECOS:\n" + pformat(solution))
416 # Retrieve primals.
417 primals = {}
418 if self.ext.options.primals is not False:
419 for variable in self.ext.variables.values():
420 offset = self._ecosVarOffset[variable]
421 dim = variable.dim
422 value = list(solution["x"][offset:offset + dim])
423 primals[variable] = value
425 # Retrieve duals.
426 duals = {}
427 if self.ext.options.duals is not False:
428 for constraint in self.ext.constraints.values():
429 if isinstance(constraint, DummyConstraint):
430 duals[constraint] = cvxopt.spmatrix(
431 [], [], [], constraint.size)
432 continue
434 dual = None
435 indices = self._ecosConIndices[constraint]
437 if isinstance(constraint, AffineConstraint):
438 if constraint.is_equality():
439 dual = list(-solution["y"][indices])
440 else:
441 dual = list(solution["z"][indices])
442 elif isinstance(constraint, SOCConstraint):
443 dual = solution["z"][indices]
444 dual = list(dual)
445 elif isinstance(constraint, RSOCConstraint):
446 dual = solution["z"][indices]
447 dual[1:] = -dual[1:]
448 alpha = dual[0] - dual[-1]
449 beta = dual[0] + dual[-1]
450 z = list(-2.0 * dual[1:-1])
451 dual = [alpha, beta] + z
452 elif isinstance(constraint, ExpConeConstraint):
453 zxy = solution["z"][indices]
454 dual = [zxy[1], zxy[2], zxy[0]]
455 else:
456 assert False, "Unexpected constraint type."
458 if type(dual) is list:
459 if len(dual) == 1:
460 dual = float(dual[0])
461 else:
462 dual = cvxopt.matrix(dual, constraint.size)
464 duals[constraint] = dual
466 # Retrieve objective value.
467 p = solution["info"]["pcost"]
468 d = solution["info"]["dcost"]
470 if p is not None and d is not None:
471 value = 0.5 * (p + d)
472 elif p is not None:
473 value = p
474 elif d is not None:
475 value = d
476 else:
477 value = None
479 if value is not None:
480 # Add back the constant part.
481 value += self.ext.no.function._constant_coef[0]
483 # Flip back the sign for maximization.
484 if self.ext.no.direction == "max":
485 value = -value
487 # Retrieve solution status.
488 statusCode = solution["info"]["exitFlag"]
489 if statusCode == 0: # optimal
490 primalStatus = SS_OPTIMAL
491 dualStatus = SS_OPTIMAL
492 problemStatus = PS_FEASIBLE
493 elif statusCode == 10: # inaccurate/optimal
494 primalStatus = SS_OPTIMAL
495 dualStatus = SS_OPTIMAL
496 problemStatus = PS_FEASIBLE
497 elif statusCode == 1: # primal infeasible
498 primalStatus = SS_INFEASIBLE
499 dualStatus = SS_UNKNOWN
500 problemStatus = PS_INFEASIBLE
501 elif statusCode == 11: # inaccurate/primal infeasible
502 primalStatus = SS_INFEASIBLE
503 dualStatus = SS_UNKNOWN
504 problemStatus = PS_INFEASIBLE
505 elif statusCode == 2: # dual infeasible
506 primalStatus = SS_UNKNOWN
507 dualStatus = SS_INFEASIBLE
508 problemStatus = PS_UNBOUNDED
509 elif statusCode == 12: # inaccurate/dual infeasible
510 primalStatus = SS_UNKNOWN
511 dualStatus = SS_INFEASIBLE
512 problemStatus = PS_UNBOUNDED
513 elif statusCode == -1: # iteration limit reached
514 primalStatus = SS_PREMATURE
515 dualStatus = SS_PREMATURE
516 problemStatus = PS_UNKNOWN
517 elif statusCode == -2: # search direction unreliable
518 primalStatus = SS_UNKNOWN
519 dualStatus = SS_UNKNOWN
520 problemStatus = PS_UNSTABLE
521 elif statusCode == -3: # numerically problematic
522 primalStatus = SS_UNKNOWN
523 dualStatus = SS_UNKNOWN
524 problemStatus = PS_UNSTABLE
525 elif statusCode == -4: # interrupted
526 primalStatus = SS_PREMATURE
527 dualStatus = SS_PREMATURE
528 problemStatus = PS_UNKNOWN
529 elif statusCode == -7: # solver error
530 primalStatus = SS_FAILURE
531 dualStatus = SS_FAILURE
532 problemStatus = PS_UNKNOWN
533 else:
534 primalStatus = SS_UNKNOWN
535 dualStatus = SS_UNKNOWN
536 problemStatus = PS_UNKNOWN
538 return self._make_solution(value, primals, duals, primalStatus,
539 dualStatus, problemStatus, {"ecos_sol": solution})
542# --------------------------------------
543__all__ = api_end(_API_START, globals())