Coverage for picos/solvers/solver_ecos.py : 85.35%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# ------------------------------------------------------------------------------
2# Copyright (C) 2018-2019 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
23from ..apidoc import api_end, api_start
24from ..constraints import (AffineConstraint, DummyConstraint,
25 ExpConeConstraint, RSOCConstraint, SOCConstraint)
26from ..expressions import AffineExpression, BinaryVariable, IntegerVariable
27from ..modeling.footprint import Specification
28from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED,
29 PS_UNKNOWN, PS_UNSTABLE, SS_FAILURE,
30 SS_INFEASIBLE, SS_OPTIMAL, SS_PREMATURE,
31 SS_UNKNOWN)
32from .solver import Solver
34_API_START = api_start(globals())
35# -------------------------------
38class ECOSSolver(Solver):
39 """Interface to the ECOS solver via its official Python interface."""
41 SUPPORTED = Specification(
42 objectives=[
43 AffineExpression],
44 constraints=[
45 DummyConstraint,
46 AffineConstraint,
47 SOCConstraint,
48 RSOCConstraint,
49 ExpConeConstraint])
51 @classmethod
52 def supports(cls, footprint, explain=False):
53 """Implement :meth:`~.solver.Solver.supports`."""
54 result = Solver.supports(footprint, explain)
55 if not result or (explain and not result[0]):
56 return result
58 if footprint not in cls.SUPPORTED:
59 if explain:
60 return False, cls.SUPPORTED.mismatch_reason(footprint)
61 else:
62 return False
64 return (True, None) if explain else True
66 @classmethod
67 def default_penalty(cls):
68 """Implement :meth:`~.solver.Solver.default_penalty`."""
69 return 1.0 # Stable free/open source solver.
71 @classmethod
72 def test_availability(cls):
73 """Implement :meth:`~.solver.Solver.test_availability`."""
74 cls.check_import("ecos")
76 @classmethod
77 def names(cls):
78 """Implement :meth:`~.solver.Solver.names`."""
79 return "ecos", "ECOS", "Embedded Conic Solver"
81 @classmethod
82 def is_free(cls):
83 """Implement :meth:`~.solver.Solver.is_free`."""
84 return True
86 def __init__(self, problem):
87 """Initialize an ECOS solver interface.
89 :param ~picos.Problem problem: The problem to be solved.
90 """
91 super(ECOSSolver, self).__init__(problem)
93 self._numVars = 0
94 """Total number of scalar variables passed to ECOS."""
96 self._ecosVarOffset = {}
97 """Maps a PICOS variable to its offset in the ECOS constraint matrix."""
99 self._ecosConIndices = dict()
100 """Maps a PICOS constraints to a range of ECOS constraint indices.
102 Note that equality constraints use a different index space."""
104 self._ecosModuleCache = None
105 """A cached reference to the ECOS module."""
107 def reset_problem(self):
108 """Implement :meth:`~.solver.Solver.reset_problem`."""
109 self.int = None
111 self._numVars = 0
112 self._ecosVarOffset.clear()
113 self._ecosConIndices.clear()
115 @property
116 def ecos(self):
117 """Return the ECOS core module ecos.py.
119 The module is obtained by ``import ecos`` up to ECOS 2.0.6 and by
120 ``import ecos.ecos`` starting with ECOS 2.0.7.
121 """
122 if self._ecosModuleCache is None:
123 import ecos
124 if hasattr(ecos, "ecos"):
125 self._ecosModuleCache = ecos.ecos
126 else:
127 self._ecosModuleCache = ecos
128 return self._ecosModuleCache
130 @property
131 def array(self):
132 """ECOS' array type."""
133 return self.ecos.np.array
135 @property
136 def matrix(self):
137 """ECOS' matrix type."""
138 return self.ecos.sparse.csc_matrix
140 def zeros(self, shape):
141 """Create a zero array or a zero matrix, depending on ``shape``."""
142 if isinstance(shape, int) or len(shape) == 1:
143 return self.ecos.np.zeros(shape)
144 else:
145 return self.matrix(shape)
147 def stack(self, *args):
148 """Stack vectors or matrices, the latter vertically."""
149 # In the case of matrices, stack vertically.
150 if isinstance(args[0], self.ecos.sparse.base.spmatrix):
151 for i in range(1, len(args)):
152 assert isinstance(args[i], self.ecos.sparse.base.spmatrix)
153 return self.ecos.sparse.vstack(args, format="csc")
155 # In the case of arrays, append them.
156 for i in range(len(args)):
157 assert isinstance(args[i], self.ecos.np.ndarray)
158 return self.ecos.np.hstack(args)
160 def _affine_expression_to_G_and_h(self, expression):
161 assert isinstance(expression, AffineExpression)
163 length = len(expression)
165 # Construct G.
166 I, J, V = [], [], []
167 for variable, coefficients in expression._linear_coefs.items():
168 ecosVarOffset = self._ecosVarOffset[variable]
170 if not isinstance(coefficients, cvxopt.spmatrix):
171 coefficients = cvxopt.sparse(coefficients)
173 I.extend(coefficients.I)
174 J.extend([ecosVarOffset + j for j in coefficients.J])
175 V.extend(coefficients.V)
176 G = self.matrix((V, (I, J)), (length, self._numVars))
178 # Construct h.
179 constant = expression._constant_coef
180 if not isinstance(constant, cvxopt.matrix):
181 constant = cvxopt.matrix(constant)
182 h = self.array(constant, dtype=float).flatten()
184 return G, h
186 _Gh = _affine_expression_to_G_and_h
188 def _import_variables(self):
189 offset = 0
190 for variable in self.ext.variables.values():
191 dim = variable.dim
193 # Mark integer variables.
194 if isinstance(variable, (BinaryVariable, IntegerVariable)):
195 ecosIndices = range(offset, offset + dim)
197 if isinstance(variable, BinaryVariable):
198 self.int["bool_vars_idx"].extend(ecosIndices)
199 else:
200 self.int["int_vars_idx"].extend(ecosIndices)
202 # Register the variable.
203 self._ecosVarOffset[variable] = offset
204 offset += dim
206 # Import bounds.
207 bounds = variable.bound_constraint
208 if bounds:
209 self._import_affine_constraint(bounds)
211 assert offset == self._numVars
213 def _append_equality(self, A, b):
214 offset = len(self.int["b"])
215 length = len(b)
217 self.int["A"] = self.stack(self.int["A"], A)
218 self.int["b"] = self.stack(self.int["b"], b)
220 return range(offset, offset + length)
222 def _append_inequality(self, G, h, typecode):
223 assert typecode in self.int["dims"]
225 # Make sure constraints are appened in the proper order.
226 if typecode == "l":
227 assert len(self.int["dims"]["q"]) == 0 \
228 and self.int["dims"]["e"] == 0
229 elif typecode == "q":
230 assert self.int["dims"]["e"] == 0
232 offset = len(self.int["h"])
233 length = len(h)
235 self.int["G"] = self.stack(self.int["G"], G)
236 self.int["h"] = self.stack(self.int["h"], h)
238 if typecode == "q":
239 self.int["dims"][typecode].append(length)
240 elif typecode == "e":
241 self.int["dims"][typecode] += 1
242 else:
243 self.int["dims"][typecode] += length
245 return range(offset, offset + length)
247 def _import_affine_constraint(self, constraint):
248 assert isinstance(constraint, AffineConstraint)
250 (G_smaller, h_smaller) = self._Gh(constraint.smaller)
251 (G_greater, h_greater) = self._Gh(constraint.greater)
252 G = G_smaller - G_greater
253 h = h_greater - h_smaller
255 if constraint.is_equality():
256 return self._append_equality(G, h)
257 else:
258 return self._append_inequality(G, h, "l")
260 def _import_socone_constraint(self, constraint):
261 assert isinstance(constraint, SOCConstraint)
263 (A, b) = self._Gh(constraint.ne)
264 (c, d) = self._Gh(constraint.ub)
266 return self._append_inequality(
267 self.stack(-c, -A), self.stack(d, b), "q")
269 def _import_rscone_constraint(self, constraint):
270 assert isinstance(constraint, RSOCConstraint)
272 (A, b) = self._Gh(constraint.ne)
273 (c1, d1) = self._Gh(constraint.ub1)
274 (c2, d2) = self._Gh(constraint.ub2)
276 return self._append_inequality(
277 self.stack(-(c1 + c2), -2.0*A, c2 - c1),
278 self.stack(d1 + d2, 2.0*b, d1 - d2), "q")
280 def _import_expcone_constraint(self, constraint):
281 assert isinstance(constraint, ExpConeConstraint)
283 (Gx, hx) = self._Gh(constraint.x)
284 (Gy, hy) = self._Gh(constraint.y)
285 (Gz, hz) = self._Gh(constraint.z)
287 # ECOS' exponential cone is cl{(x,y,z) | exp(x/z) <= y/z, z > 0},
288 # PICOS' is cl{(x,y,z) | x >= y*exp(z/y), y > 0}. Note that given y > 0
289 # it is x >= y*exp(z/y) if and only if exp(z/y) <= x/y. Therefor we can
290 # transform from our coordinates to theirs with the mapping
291 # (x, y, z) ↦ (z, x, y). Further, G and h with G = (Gx, Gy, Gz) and
292 # h = (hx, hy, hz) are such that G*X + h = (x, y, z) where X is the
293 # row-vectorization of all variables. ECOS however expects G and h such
294 # that h - G*X is constrained to be in the exponential cone.
295 return self._append_inequality(
296 -self.stack(Gz, Gx, Gy), self.stack(hz, hx, hy), "e")
298 def _import_objective(self):
299 direction, objective = self.ext.no
301 # ECOS only supports minimization; flip the sign for maximization.
302 sign = 1.0 if direction == "min" else -1.0
304 # Import coefficients.
305 for variable, coefficient in objective._linear_coefs.items():
306 for localIndex in range(variable.dim):
307 index = self._ecosVarOffset[variable] + localIndex
308 self.int["c"][index] = sign * coefficient[localIndex]
310 def _import_problem(self):
311 self._numVars = sum(var.dim for var in self.ext.variables.values())
313 # ECOS' internal problem representation is stateless; a number of
314 # vectors and matrices is supplied each time a search is started.
315 # These vectors and matrices are thus stored in self.int.
316 self.int = {
317 # Objective function coefficients.
318 "c": self.zeros(self._numVars),
320 # Linear equality left hand side.
321 "A": self.matrix((0, self._numVars)),
323 # Linear equality right hand side.
324 "b": self.array([]),
326 # Conic inequality left hand side.
327 "G": self.matrix((0, self._numVars)),
329 # Conic inequality right hand side.
330 "h": self.array([]),
332 # Cone definition: Linear, second order, exponential dimensions.
333 "dims": {"l": 0, "q": [], "e": 0},
335 # Boolean variable indices.
336 "bool_vars_idx": [],
338 # Integer variable indices.
339 "int_vars_idx": []
340 }
342 # NOTE: Constraints need to be append ordered by type.
343 # TODO: Add fast constraints-by-type iterators to Problem.
345 # Import variables with their bounds as affine constraints.
346 self._import_variables()
348 # Import affine constraints.
349 for constraint in self.ext.constraints.values():
350 if isinstance(constraint, AffineConstraint):
351 self._ecosConIndices[constraint] = \
352 self._import_affine_constraint(constraint)
354 # Import second order cone constraints.
355 for constraint in self.ext.constraints.values():
356 if isinstance(constraint, SOCConstraint):
357 self._ecosConIndices[constraint] = \
358 self._import_socone_constraint(constraint)
360 # Import rotated second order cone constraints.
361 for constraint in self.ext.constraints.values():
362 if isinstance(constraint, RSOCConstraint):
363 self._ecosConIndices[constraint] = \
364 self._import_rscone_constraint(constraint)
366 # Import exponential cone constraints.
367 for constraint in self.ext.constraints.values():
368 if isinstance(constraint, ExpConeConstraint):
369 self._ecosConIndices[constraint] = \
370 self._import_expcone_constraint(constraint)
372 # Make sure that no unsupported constraints are present.
373 for constraint in self.ext.constraints.values():
374 assert isinstance(constraint, (DummyConstraint, AffineConstraint,
375 SOCConstraint, RSOCConstraint, ExpConeConstraint))
377 # Set objective.
378 self._import_objective()
380 def _update_problem(self):
381 raise NotImplementedError
383 def _solve(self):
384 ecosOptions = {}
386 # verbosity
387 beVerbose = self.verbosity() >= 1
388 ecosOptions["verbose"] = beVerbose
389 ecosOptions["mi_verbose"] = beVerbose
391 # rel_prim_fsb_tol, rel_dual_fsb_tol
392 feasibilityTols = [tol for tol in (self.ext.options.rel_prim_fsb_tol,
393 self.ext.options.rel_dual_fsb_tol) if tol is not None]
394 if feasibilityTols:
395 ecosOptions["feastol"] = min(feasibilityTols)
397 # abs_ipm_opt_tol
398 if self.ext.options.abs_ipm_opt_tol is not None:
399 ecosOptions["abstol"] = self.ext.options.abs_ipm_opt_tol
401 # rel_ipm_opt_tol
402 if self.ext.options.rel_ipm_opt_tol is not None:
403 ecosOptions["reltol"] = self.ext.options.rel_ipm_opt_tol
405 # abs_bnb_opt_tol
406 if self.ext.options.abs_bnb_opt_tol is not None:
407 ecosOptions["mi_abs_eps"] = self.ext.options.abs_bnb_opt_tol
409 # rel_bnb_opt_tol
410 if self.ext.options.rel_bnb_opt_tol is not None:
411 ecosOptions["mi_rel_eps"] = self.ext.options.rel_bnb_opt_tol
413 # integrality_tol
414 # HACK: ECOS_BB uses ECOS' "tolerance for feasibility condition for
415 # inaccurate solution" as the integrality tolerance.
416 if self.ext.options.integrality_tol is not None:
417 ecosOptions["feastol_inacc"] = self.ext.options.integrality_tol
419 # max_iterations
420 if self.ext.options.max_iterations is not None:
421 ecosOptions["max_iters"] = self.ext.options.max_iterations
422 ecosOptions["mi_max_iters"] = self.ext.options.max_iterations
424 # Handle unsupported options.
425 self._handle_unsupported_options(
426 "lp_root_method", "lp_node_method", "timelimit", "treememory",
427 "max_fsb_nodes", "hotstart")
429 # Assemble the solver input.
430 arguments = {}
431 arguments.update(self.int)
432 arguments.update(ecosOptions)
434 # Debug print the solver input.
435 if self._debug():
436 from pprint import pformat
437 self._debug("Passing to ECOS:\n" + pformat(arguments))
439 # Attempt to solve the problem.
440 with self._header(), self._stopwatch():
441 solution = self.ecos.solve(**arguments)
443 # Debug print the solver output.
444 if self._debug():
445 from pprint import pformat
446 self._debug("Recevied from ECOS:\n" + pformat(solution))
448 # Retrieve primals.
449 primals = {}
450 if self.ext.options.primals is not False:
451 for variable in self.ext.variables.values():
452 offset = self._ecosVarOffset[variable]
453 dim = variable.dim
454 value = list(solution["x"][offset:offset + dim])
455 primals[variable] = value
457 # Retrieve duals.
458 duals = {}
459 if self.ext.options.duals is not False:
460 for constraint in self.ext.constraints.values():
461 if isinstance(constraint, DummyConstraint):
462 duals[constraint] = cvxopt.spmatrix(
463 [], [], [], constraint.size)
464 continue
466 dual = None
467 indices = self._ecosConIndices[constraint]
469 if isinstance(constraint, AffineConstraint):
470 if constraint.is_equality():
471 dual = list(-solution["y"][indices])
472 else:
473 dual = list(solution["z"][indices])
474 elif isinstance(constraint, SOCConstraint):
475 dual = solution["z"][indices]
476 dual = list(dual)
477 elif isinstance(constraint, RSOCConstraint):
478 dual = solution["z"][indices]
479 dual[1:] = -dual[1:]
480 alpha = dual[0] - dual[-1]
481 beta = dual[0] + dual[-1]
482 z = list(-2.0 * dual[1:-1])
483 dual = [alpha, beta] + z
484 elif isinstance(constraint, ExpConeConstraint):
485 zxy = solution["z"][indices]
486 dual = [zxy[1], zxy[2], zxy[0]]
487 else:
488 assert False, "Unexpected constraint type."
490 if type(dual) is list:
491 if len(dual) == 1:
492 dual = float(dual[0])
493 else:
494 dual = cvxopt.matrix(dual, constraint.size)
496 duals[constraint] = dual
498 # Retrieve objective value.
499 p = solution["info"]["pcost"]
500 d = solution["info"]["dcost"]
502 if p is not None and d is not None:
503 value = 0.5 * (p + d)
504 elif p is not None:
505 value = p
506 elif d is not None:
507 value = d
508 else:
509 value = None
511 if value is not None:
512 # Add back the constant part.
513 value += self.ext.no.function._constant_coef[0]
515 # Flip back the sign for maximization.
516 if self.ext.no.direction == "max":
517 value = -value
519 # Retrieve solution status.
520 statusCode = solution["info"]["exitFlag"]
521 if statusCode == 0: # optimal
522 primalStatus = SS_OPTIMAL
523 dualStatus = SS_OPTIMAL
524 problemStatus = PS_FEASIBLE
525 elif statusCode == 10: # inaccurate/optimal
526 primalStatus = SS_OPTIMAL
527 dualStatus = SS_OPTIMAL
528 problemStatus = PS_FEASIBLE
529 elif statusCode == 1: # primal infeasible
530 primalStatus = SS_INFEASIBLE
531 dualStatus = SS_UNKNOWN
532 problemStatus = PS_INFEASIBLE
533 elif statusCode == 11: # inaccurate/primal infeasible
534 primalStatus = SS_INFEASIBLE
535 dualStatus = SS_UNKNOWN
536 problemStatus = PS_INFEASIBLE
537 elif statusCode == 2: # dual infeasible
538 primalStatus = SS_UNKNOWN
539 dualStatus = SS_INFEASIBLE
540 problemStatus = PS_UNBOUNDED
541 elif statusCode == 12: # inaccurate/dual infeasible
542 primalStatus = SS_UNKNOWN
543 dualStatus = SS_INFEASIBLE
544 problemStatus = PS_UNBOUNDED
545 elif statusCode == -1: # iteration limit reached
546 primalStatus = SS_PREMATURE
547 dualStatus = SS_PREMATURE
548 problemStatus = PS_UNKNOWN
549 elif statusCode == -2: # search direction unreliable
550 primalStatus = SS_UNKNOWN
551 dualStatus = SS_UNKNOWN
552 problemStatus = PS_UNSTABLE
553 elif statusCode == -3: # numerically problematic
554 primalStatus = SS_UNKNOWN
555 dualStatus = SS_UNKNOWN
556 problemStatus = PS_UNSTABLE
557 elif statusCode == -4: # interrupted
558 primalStatus = SS_PREMATURE
559 dualStatus = SS_PREMATURE
560 problemStatus = PS_UNKNOWN
561 elif statusCode == -7: # solver error
562 primalStatus = SS_FAILURE
563 dualStatus = SS_FAILURE
564 problemStatus = PS_UNKNOWN
565 else:
566 primalStatus = SS_UNKNOWN
567 dualStatus = SS_UNKNOWN
568 problemStatus = PS_UNKNOWN
570 return self._make_solution(value, primals, duals, primalStatus,
571 dualStatus, problemStatus, {"ecos_sol": solution})
574# --------------------------------------
575__all__ = api_end(_API_START, globals())