Coverage for picos/expressions/uncertain/pert_conic.py: 88.89%
180 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) 2020 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"""Implements :class:`ConicPerturbationSet`."""
21from collections import namedtuple
23import cvxopt
25from ... import glyphs, settings
26from ...apidoc import api_end, api_start
27from ...caching import cached_property
28from ...constraints import SOCConstraint
29from ..cone_product import ProductCone
30from ..cone_soc import SecondOrderCone
31from ..data import cvxopt_hcat, cvxopt_vcat
32from ..exp_affine import Constant
33from ..expression import NotValued
34from ..set_ball import Ball
35from ..variables import RealVariable
36from ..vectorizations import FullVectorization
37from .perturbation import Perturbation, PerturbationUniverse
38from .uexp_affine import UncertainAffineExpression
39from .uexpression import IntractableWorstCase
41_API_START = api_start(globals())
42# -------------------------------
45class ConicPerturbationSet(PerturbationUniverse):
46 r"""A conic description of a :class:`~.perturbation.Perturbation`.
48 :Definition:
50 An instance :math:`\Theta` of this class defines a perturbation parameter
52 .. math::
54 \theta \in \Theta = \{t \in \mathbb{R}^{m \times n}
55 \mid \exists u \in \mathbb{R}^l
56 : A\operatorname{vec}(t) + Bu + c \in K\}
58 where :math:`m, n \in \mathbb{Z}_{\geq 1}`,
59 :math:`l \in \mathbb{Z}_{\geq 0}`, :math:`A \in \mathbb{R}^{k \times mn}`,
60 :math:`B \in \mathbb{R}^{k \times l}`, :math:`c \in \mathbb{R}^k` and
61 :math:`K \subseteq \mathbb{R}^k` is a (product) cone for some
62 :math:`k \in \mathbb{Z}_{\geq 1}`.
64 :Usage:
66 Obtaining :math:`\theta` is done in a number of steps:
68 1. Create an instance of this class (see :meth:`__init__`).
69 2. Access :attr:`element` to obtain a regular, fresh
70 :class:`~.variables.RealVariable` representing :math:`t`.
71 3. Define :math:`\Theta` through any number of regular PICOS constraints
72 that depend only on :math:`t` and that have a conic representation by
73 passing the constraints to :meth:`bound`.
74 4. Call :meth:`compile` to obtain a handle to the
75 :class:`~.perturbation.Perturbation` :math:`\theta`.
76 5. You can now use :math:`\theta` to build instances of
77 :class:`~.uexp_affine.UncertainAffineExpression` and derived constraint
78 types.
80 It is best practice to assign :math:`t` to a Python variable and overwrite
81 that variable with :math:`\theta` on compilation.
83 Alternatively, you can obtain a compiled :math:`\Theta` from the factory
84 method :meth:`from_constraints` and access :math:`\theta` via
85 :attr:`parameter`.
87 :Example:
89 >>> from picos import Constant, Norm, RealVariable
90 >>> from picos.uncertain import ConicPerturbationSet
91 >>> S = ConicPerturbationSet("P", (4, 4))
92 >>> P = S.element; P # Obtain a temporary parameter to describe S with.
93 <4×4 Real Variable: P>
94 >>> S.bound(Norm(P, float("inf")) <= 1) # Confine each element to [-1,1].
95 >>> S.bound(Norm(P, 1) <= 4); S # Allow a perturbation budget of 4.
96 <4×4 Conic Perturbation Set: {P : ‖P‖_max ≤ 1 ∧ ‖P‖_sum ≤ 4}>
97 >>> P = S.compile(); P # Compile the set and obtain the actual parameter.
98 <4×4 Perturbation: P>
99 >>> A = Constant("A", range(16), (4, 4))
100 >>> U = A + P; U # Define an uncertain data matrix.
101 <4×4 Uncertain Affine Expression: A + P>
102 >>> x = RealVariable("x", 4)
103 >>> U*x # Define an uncertain affine expression.
104 <4×1 Uncertain Affine Expression: (A + P)·x>
105 """
107 def __init__(self, parameter_name, shape=(1, 1)):
108 """Create a :class:`ConicPerturbationSet`.
110 :param str parameter_name:
111 Name of the parameter that lives in the set.
113 :param shape:
114 The shape of a vector or matrix perturbation.
115 :type shape:
116 int or tuple or list
117 """
118 from ...modeling import Problem
120 self._compiled = False
121 self._element = RealVariable(parameter_name, shape)
122 self._parameter = Perturbation(self, parameter_name, shape)
123 self._bounds = Problem()
125 Subtype = namedtuple("Subtype", (
126 "param_dim", "cone_type", "dual_cone_type", "has_B"))
128 def _subtype(self):
129 return self.Subtype(
130 self._parameter.dim,
131 self.K.type,
132 self.K.dual_cone.type,
133 self.B is not None)
135 @classmethod
136 def from_constraints(cls, parameter_name, *constraints):
137 """Create a :class:`ConicPerturbationSet` from constraints.
139 The constraints must concern a single regular decision variable that
140 plays the role of the :attr:`element` :math:`t`. This variable is not
141 stored or modified and can be reused in a different context.
143 :param str parameter_name:
144 Name of the parameter that lives in the set.
146 :param constraints:
147 A parameter sequence of constraints that concern a single regular
148 decision variable whose internal vectorization is trivial (its
149 dimension must match the product over its shape) and that have a
150 conic representation.
152 :raises ValueError:
153 If the constraints do not all concern the same single variable.
155 :raises TypeError:
156 If the variable uses a nontrivial vectorization format or if the
157 constraints do not all have a conic representation.
159 :Example:
161 >>> from picos.expressions.uncertain import ConicPerturbationSet
162 >>> from picos import RealVariable
163 >>> x = RealVariable("x", 4)
164 >>> T = ConicPerturbationSet.from_constraints("t", abs(x) <= 2, x >= 0)
165 >>> print(T)
166 {t : ‖t‖ ≤ 2 ∧ t ≥ 0}
167 >>> print(repr(T.parameter))
168 <4×1 Perturbation: t>
169 """
170 T, t, seen_variable = None, None, None
172 for constraint in constraints:
173 if len(constraint.variables) != 1:
174 raise ValueError("The constraint {} does not concern exactly "
175 "one variable.".format(constraint))
177 variable = next(iter(constraint.variables))
179 if not isinstance(variable._vec, FullVectorization):
180 raise TypeError("The variable {} cannot be used to construct a "
181 "{} from constraints because it uses a nontrivial "
182 "vectorization format.".format(variable, cls.name))
184 if not seen_variable:
185 seen_variable = variable
186 T = cls(parameter_name, variable.shape)
187 t = T.element
188 elif variable is not seen_variable:
189 raise ValueError("The constraints do not concern the same "
190 "single variable (found {} and {})."
191 .format(seen_variable.name, variable.name))
193 T.bound(constraint.replace_mutables({variable: t}))
195 T.compile()
196 return T
198 def __str__(self):
199 return glyphs.set(glyphs.sep(self._parameter.name,
200 glyphs.and_("", "").join(str(con)
201 for con in self._bounds.constraints.values())))
203 @classmethod
204 def _get_type_string_base(cls):
205 return "Conic Perturbation Set"
207 def __repr__(self):
208 return glyphs.repr2("{} {}".format(glyphs.shape(self._element.shape),
209 self._get_type_string_base()), self.__str__())
211 def _forbid_compiled(self):
212 if self._compiled:
213 raise RuntimeError(
214 "The perturbation set has already been compiled.")
216 def _require_compiled(self):
217 if not self._compiled:
218 raise RuntimeError(
219 "The perturbation set has not yet been compiled.")
221 @property
222 def element(self):
223 r"""The perturbation element :math:`t` describing the set.
225 This is a regular :class:`~.variables.RealVariable` that you can use to
226 create constraints to pass to :meth:`bound`. You can then obtain the
227 "actual" perturbation parameter :math:`\theta` to use in expressions
228 alongside your decision variaiables using :meth:`compile`.
230 .. warning::
232 If you use this object instead of :attr:`parameter` to define a
233 decision problem then it will act as a regular decision variable,
234 which is probably not what you want.
236 :raises RuntimeError:
237 If the set was already compiled.
238 """
239 self._forbid_compiled()
240 return self._element
242 def bound(self, constraint):
243 r"""Add a constraint that bounds :math:`t`.
245 The constraints do not need to be conic but they need to have a *conic
246 representation*, which may involve any number of auxiliary variables.
247 For instance, given a constant *uncertainty budget* :math:`b`, you may
248 add the bound :math:`\lVert t \rVert_1 \leq b` (via
249 ``picos.Norm(t, 1) <= b``) which can be represented in conic form as
251 .. math::
253 &\exists v \in \mathbb{R}^{\operatorname{dim}(t)}
254 : -t \leq v \land t \leq v \land \mathbf{1}^Tv \leq b \\
255 \Longleftrightarrow~
256 &\exists v \in \mathbb{R}^{\operatorname{dim}(t)} :
257 \begin{pmatrix}
258 v + t \\
259 v - t \\
260 b - \mathbf{1}^Tv
261 \end{pmatrix} \in \mathbb{R}_{\geq 0}^{2\operatorname{dim}(t) + 1}.
263 The auxiliary variable :math:`v` then plays the role of (a slice of)
264 :math:`u` in the formal definition of :math:`\Theta`.
266 When you are done adding bounds, you can obtain :math:`\theta` using
267 :meth:`compile`.
269 :raises RuntimeError:
270 If the set was already compiled.
271 """
272 self._forbid_compiled()
274 vars = constraint.variables
276 if len(vars) != 1 or next(iter(vars)) != self._element:
277 raise ValueError(
278 "The constraint {} does not bound (only) the perturbation "
279 "element {}.".format(constraint, self._element.string))
281 self._bounds.add_constraint(constraint)
283 def compile(self, validate_feasibility=False):
284 r"""Compile the set and return :math:`\theta`.
286 Internally, this computes the matrices :math:`A` and :math:`B`, the
287 vector :math:`c` and the (product) cone :math:`K`.
289 :param bool validate_feasibility:
290 Whether to solve the feasibility problem associated with the bounds
291 on :math:`t` to verify that :math:`\Theta` is nonempty.
293 :returns:
294 An instance of :class:`~.perturbation.Perturbation`.
296 :raises RuntimeError:
297 If the set was already compiled.
299 :raises TypeError:
300 If the bound constraints could not be put into conic form.
302 :raises ValueError:
303 If :math:`\Theta` could not be verified to be nonempty (needs
304 ``validate_feasibility=True``).
305 """
306 from ...constraints import DummyConstraint
307 from ...modeling import NoStrategyFound
309 self._forbid_compiled()
311 # If no bounds are given, add a DummyConstraint to mark t free.
312 if not self._bounds.constraints:
313 self._bounds.add_constraint(DummyConstraint(self._element))
315 # Transform all bounds into conic form.
316 try:
317 conic_bounds = self._bounds.conic_form
318 except NoStrategyFound as error:
319 raise TypeError("Could not find a conic representation for all of "
320 "the bounds on {}.".format(self._element.string)) from error
322 # Validate bound feasibility if requested.
323 if validate_feasibility:
324 from ...modeling.solution import SS_OPTIMAL, SS_FEASIBLE
326 solution = conic_bounds.solve(
327 primals=None, apply_solution=False, **settings.INTERNAL_OPTIONS)
329 status = solution.primalStatus
331 if status not in (SS_OPTIMAL, SS_FEASIBLE):
332 raise ValueError("Could not verify that the bounds on {} are "
333 "feasible: The solver {} reports a primal solution state of"
334 " {} for the associated feasibility problem.".format(
335 self._element.string, solution.solver, status))
337 # Form a virtual variable u from the auxiliary variables.
338 u = [var for var in conic_bounds.variables.values()
339 if var is not self._element]
341 # Convert auxiliary variable bounds to additional affine constraints.
342 conic_bounds.add_list_of_constraints([
343 var.bound_constraint for var in u if var.bound_constraint])
345 # Reformulate bounds with respect to a single product cone K.
346 A, B, c, K = [], [], [], []
347 for constraint in conic_bounds.constraints.values():
348 member, cone = constraint.conic_membership_form
350 if self._element in member._linear_coefs:
351 A.append(member._linear_coefs[self._element])
352 else:
353 A.append(cvxopt.spmatrix(
354 [], [], [], size=(cone.dim, self._element.dim)))
356 if u:
357 B.append(cvxopt_hcat([
358 member._linear_coefs[var] if var in member._linear_coefs
359 else cvxopt.spmatrix([], [], [], size=(cone.dim, var.dim))
360 for var in u]))
362 c.append(member._constant_coef)
363 K.append(cone)
365 self._A = Constant("A", cvxopt_vcat(A))
366 self._B = Constant("B", cvxopt_vcat(B)) if u else None
367 self._c = Constant("c", cvxopt_vcat(c))
368 self._K = ProductCone(*K)
370 # Store the bounds in conic form to speed up worst_case.
371 self._conic_bounds = conic_bounds
373 self._compiled = True
374 return self._parameter
376 @property
377 def distributional(self):
378 """Implement for :class:`~.perturbation.PerturbationUniverse`."""
379 return False
381 @property
382 def parameter(self):
383 r"""The perturbation parameter :math:`\theta` living in the set.
385 This is the object returned by :meth:`compile`.
387 :raises RuntimeError:
388 If the set has not been compiled.
389 """
390 self._require_compiled()
391 return self._parameter
393 @property
394 def A(self):
395 r"""The compiled matrix :math:`A`.
397 :raises RuntimeError:
398 If the set has not been compiled.
399 """
400 self._require_compiled()
401 return self._A
403 @property
404 def B(self):
405 r"""The compiled matrix :math:`B` or :obj:`None` if :math:`l = 0`.
407 :raises RuntimeError:
408 If the set has not been compiled.
409 """
410 self._require_compiled()
411 return self._B
413 @property
414 def c(self):
415 r"""The compiled vector :math:`c`.
417 :raises RuntimeError:
418 If the set has not been compiled.
419 """
420 self._require_compiled()
421 return self._c
423 @property
424 def K(self):
425 r"""The compiled (product) cone :math:`K`.
427 :raises RuntimeError:
428 If the set has not been compiled.
429 """
430 self._require_compiled()
431 return self._K
433 def worst_case(self, scalar, direction):
434 """Implement for :class:`~.perturbation.PerturbationUniverse`."""
435 from ...modeling import SolutionFailure
437 self._require_compiled()
438 self._check_worst_case_argument_scalar(scalar)
439 self._check_worst_case_argument_direction(direction)
441 p = self._parameter
442 P = self._conic_bounds.copy()
443 x = P.variables[p.name]
444 f = scalar.replace_mutables({p: x})
446 self._check_worst_case_f_and_x(f, x)
448 if (direction == "min" and not f.convex) \
449 or (direction == "max" and not f.concave):
450 raise IntractableWorstCase("PICOS refuses to compute {}({}) for {} "
451 "as this is a nonconvex problem.".format(direction, f.string,
452 glyphs.element(p.name, self)))
454 P.set_objective(direction, f if direction != "find" else None)
456 try:
457 P.solve(**settings.INTERNAL_OPTIONS)
458 return f.safe_value, x.safe_value
459 except (SolutionFailure, NotValued) as error:
460 raise RuntimeError("Failed to compute {}({}) for {}.".format(
461 direction, f.string, glyphs.element(x.string, self))) from error
463 @cached_property
464 def unit_ball_form(self):
465 """A recipe to repose from ellipsoidal to unit norm ball uncertainty.
467 If the set is :attr:`ellipsoidal`, then this is a pair ``(U, M)`` where
468 ``U`` is a :class:`~.pert_conic.UnitBallPerturbationSet` and ``M`` is a
469 dictionary mapping the old :attr:`parameter` to an affine expression of
470 the new parameter that can represent the old parameter in an expression
471 (see :meth:`~.expression.Expression.replace_mutables`).
472 The mapping ``M`` is empty if and only if the perturbation set is
473 already an instance of :class:`~.pert_conic.UnitBallPerturbationSet`.
475 If the uncertainty set is not ellipsoidal, then this is :obj:`None`.
477 See also :attr:`SOCConstraint.unit_ball_form
478 <.con_soc.SOCConstraint.unit_ball_form>`.
480 :Example:
482 >>> from picos import Problem, RealVariable, sum
483 >>> from picos.uncertain import ConicPerturbationSet
484 >>> # Create a conic perturbation set and a refinement recipe.
485 >>> T = ConicPerturbationSet("t", (2, 2))
486 >>> T.bound(abs(([[1, 2], [3, 4]] ^ T.element) + 1) <= 10)
487 >>> t = T.compile()
488 >>> U, mapping = T.unit_ball_form
489 >>> print(U)
490 {t' : ‖t'‖ ≤ 1}
491 >>> print(mapping)
492 {<2×2 Perturbation: t>: <2×2 Uncertain Affine Expression: t(t')>}
493 >>> # Define and solve a conically uncertain LP.
494 >>> X = RealVariable("X", (2, 2))
495 >>> P = Problem()
496 >>> P.set_objective("max", sum(X))
497 >>> _ = P.add_constraint(X + 2*t <= 10)
498 >>> print(repr(P.parameters["t"].universe))
499 <2×2 Conic Perturbation Set: {t : ‖[2×2]⊙t + [1]‖ ≤ 10}>
500 >>> _ = P.solve(solver="cvxopt")
501 >>> print(X)
502 [-8.00e+00 1.00e+00]
503 [ 4.00e+00 5.50e+00]
504 >>> # Refine the problem to a unit ball uncertain LP.
505 >>> Q = Problem()
506 >>> Q.set_objective("max", sum(X))
507 >>> _ = Q.add_constraint(X + 2*mapping[t] <= 10)
508 >>> print(repr(Q.parameters["t'"].universe))
509 <2×2 Unit Ball Perturbation Set: {t' : ‖t'‖ ≤ 1}>
510 >>> _ = Q.solve(solver="cvxopt")
511 >>> print(X)
512 [-8.00e+00 1.00e+00]
513 [ 4.00e+00 5.50e+00]
514 """
515 self._require_compiled()
517 if self._B is not None:
518 return None
520 K = self._K.refined
522 if not isinstance(K, SecondOrderCone):
523 return None
525 C = self._A*self._element.vec + self._c << K
527 assert isinstance(C, SOCConstraint)
529 try:
530 X, aff_y, y, _ = C.unit_ball_form
531 except ValueError:
532 return None
534 assert X is self._element
536 U = UnitBallPerturbationSet("{}'".format(X.name), X.shape)
537 u = U.parameter
538 replacement = UncertainAffineExpression("{}({})".format(X.name, u.name),
539 X.shape, {(u,): aff_y._linear_coefs[y], (): aff_y._constant_coef})
541 return U, {self._parameter: replacement}
543 @property
544 def ellipsoidal(self):
545 """Whether the perturbation set is an ellipsoid.
547 If this is true, then a :attr:`unit_ball_form` is available.
548 """
549 return bool(self.unit_ball_form)
552class UnitBallPerturbationSet(ConicPerturbationSet):
553 r"""Represents perturbation in an Euclidean or Frobenius unit norm ball.
555 This is a :class:`~.pert_conic.ConicPerturbationSet` with fixed form
557 .. math::
559 \{t \in \mathbb{R}^{m \times n} \mid \lVert t \rVert_F \leq 1\}.
561 After initialization, you can obtain the parameter using
562 :attr:`~.pert_conic.ConicPerturbationSet.parameter`.
563 """
565 def __init__(self, parameter_name, shape=(1, 1)):
566 """See :meth:`ConicPerturbationSet.__init__`."""
567 ConicPerturbationSet.__init__(self, parameter_name, shape)
568 self.bound(self._element << Ball())
569 self.compile()
571 @classmethod
572 def _get_type_string_base(cls):
573 return "Unit Ball Perturbation Set"
575 @property
576 def unit_ball_form(self):
577 """Overwrite :attr:`ConicPerturbationSet.unit_ball_form`."""
578 return self, {}
581# --------------------------------------
582__all__ = api_end(_API_START, globals())