Source code for picos.expressions.uncertain.pert_conic

# coding: utf-8

# ------------------------------------------------------------------------------
# Copyright (C) 2020 Maximilian Stahlberg
# This file is part of PICOS.
# PICOS is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# PICOS is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <>.
# ------------------------------------------------------------------------------

"""Implements :class:`ConicPerturbationSet`."""

from collections import namedtuple

import cvxopt

from ... import glyphs, settings
from ...apidoc import api_end, api_start
from ...caching import cached_property
from ...constraints import SOCConstraint
from ..cone_product import ProductCone
from ..cone_soc import SecondOrderCone
from import cvxopt_hcat, cvxopt_vcat
from ..exp_affine import Constant
from ..expression import NotValued
from ..set_ball import Ball
from ..variables import RealVariable
from ..vectorizations import FullVectorization
from .perturbation import Perturbation, PerturbationUniverse
from .uexp_affine import UncertainAffineExpression
from .uexpression import IntractableWorstCase

_API_START = api_start(globals())
# -------------------------------

[docs]class ConicPerturbationSet(PerturbationUniverse): r"""A conic description of a :class:`~.perturbation.Perturbation`. :Definition: An instance :math:`\Theta` of this class defines a perturbation parameter .. math:: \theta \in \Theta = \{t \in \mathbb{R}^{m \times n} \mid \exists u \in \mathbb{R}^l : A\operatorname{vec}(t) + Bu + c \in K\} where :math:`m, n \in \mathbb{Z}_{\geq 1}`, :math:`l \in \mathbb{Z}_{\geq 0}`, :math:`A \in \mathbb{R}^{k \times mn}`, :math:`B \in \mathbb{R}^{k \times l}`, :math:`c \in \mathbb{R}^k` and :math:`K \subseteq \mathbb{R}^k` is a (product) cone for some :math:`k \in \mathbb{Z}_{\geq 1}`. :Usage: Obtaining :math:`\theta` is done in a number of steps: 1. Create an instance of this class (see :meth:`__init__`). 2. Access :attr:`element` to obtain a regular, fresh :class:`~.variables.RealVariable` representing :math:`t`. 3. Define :math:`\Theta` through any number of regular PICOS constraints that depend only on :math:`t` and that have a conic representation by passing the constraints to :meth:`bound`. 4. Call :meth:`compile` to obtain a handle to the :class:`~.perturbation.Perturbation` :math:`\theta`. 5. You can now use :math:`\theta` to build instances of :class:`~.uexp_affine.UncertainAffineExpression` and derived constraint types. It is best practice to assign :math:`t` to a Python variable and overwrite that variable with :math:`\theta` on compilation. Alternatively, you can obtain a compiled :math:`\Theta` from the factory method :meth:`from_constraints` and access :math:`\theta` via :attr:`parameter`. :Example: >>> from picos import Constant, Norm, RealVariable >>> from picos.uncertain import ConicPerturbationSet >>> S = ConicPerturbationSet("P", (4, 4)) >>> P = S.element; P # Obtain a temporary parameter to describe S with. <4×4 Real Variable: P> >>> S.bound(Norm(P, float("inf")) <= 1) # Confine each element to [-1,1]. >>> S.bound(Norm(P, 1) <= 4); S # Allow a perturbation budget of 4. <4×4 Conic Perturbation Set: {P : ‖P‖_max ≤ 1 ∧ ‖P‖_sum ≤ 4}> >>> P = S.compile(); P # Compile the set and obtain the actual parameter. <4×4 Perturbation: P> >>> A = Constant("A", range(16), (4, 4)) >>> U = A + P; U # Define an uncertain data matrix. <4×4 Uncertain Affine Expression: A + P> >>> x = RealVariable("x", 4) >>> U*x # Define an uncertain affine expression. <4×1 Uncertain Affine Expression: (A + P)·x> """
[docs] def __init__(self, parameter_name, shape=(1, 1)): """Create a :class:`ConicPerturbationSet`. :param str parameter_name: Name of the parameter that lives in the set. :param shape: The shape of a vector or matrix perturbation. :type shape: int or tuple or list """ from ...modeling import Problem self._compiled = False self._element = RealVariable(parameter_name, shape) self._parameter = Perturbation(self, parameter_name, shape) self._bounds = Problem()
Subtype = namedtuple("Subtype", ( "param_dim", "cone_type", "dual_cone_type", "has_B")) def _subtype(self): return self.Subtype( self._parameter.dim, self.K.type, self.K.dual_cone.type, self.B is not None)
[docs] @classmethod def from_constraints(cls, parameter_name, *constraints): """Create a :class:`ConicPerturbationSet` from constraints. The constraints must concern a single regular decision variable that plays the role of the :attr:`element` :math:`t`. This variable is not stored or modified and can be reused in a different context. :param str parameter_name: Name of the parameter that lives in the set. :param constraints: A parameter sequence of constraints that concern a single regular decision variable whose internal vectorization is trivial (its dimension must match the product over its shape) and that have a conic representation. :raises ValueError: If the constraints do not all concern the same single variable. :raises TypeError: If the variable uses a nontrivial vectorization format or if the constraints do not all have a conic representation. :Example: >>> from picos.expressions.uncertain import ConicPerturbationSet >>> from picos import RealVariable >>> x = RealVariable("x", 4) >>> T = ConicPerturbationSet.from_constraints("t", abs(x) <= 2, x >= 0) >>> print(T) {t : ‖t‖ ≤ 2 ∧ t ≥ 0} >>> print(repr(T.parameter)) <4×1 Perturbation: t> """ T, t, seen_variable = None, None, None for constraint in constraints: if len(constraint.variables) != 1: raise ValueError("The constraint {} does not concern exactly " "one variable.".format(constraint)) variable = next(iter(constraint.variables)) if not isinstance(variable._vec, FullVectorization): raise TypeError("The variable {} cannot be used to construct a " "{} from constraints because it uses a nontrivial " "vectorization format.".format(variable, if not seen_variable: seen_variable = variable T = cls(parameter_name, variable.shape) t = T.element elif variable is not seen_variable: raise ValueError("The constraints do not concern the same " "single variable (found {} and {})." .format(, T.bound(constraint.replace_mutables({variable: t})) T.compile() return T
def __str__(self): return glyphs.set(glyphs.sep(, glyphs.and_("", "").join(str(con) for con in self._bounds.constraints.values()))) @classmethod def _get_type_string_base(cls): return "Conic Perturbation Set" def __repr__(self): return glyphs.repr2("{} {}".format(glyphs.shape(self._element.shape), self._get_type_string_base()), self.__str__()) def _forbid_compiled(self): if self._compiled: raise RuntimeError( "The perturbation set has already been compiled.") def _require_compiled(self): if not self._compiled: raise RuntimeError( "The perturbation set has not yet been compiled.") @property def element(self): r"""The perturbation element :math:`t` describing the set. This is a regular :class:`~.variables.RealVariable` that you can use to create constraints to pass to :meth:`bound`. You can then obtain the "actual" perturbation parameter :math:`\theta` to use in expressions alongside your decision variaiables using :meth:`compile`. .. warning:: If you use this object instead of :attr:`parameter` to define a decision problem then it will act as a regular decision variable, which is probably not what you want. :raises RuntimeError: If the set was already compiled. """ self._forbid_compiled() return self._element
[docs] def bound(self, constraint): r"""Add a constraint that bounds :math:`t`. The constraints do not need to be conic but they need to have a *conic representation*, which may involve any number of auxiliary variables. For instance, given a constant *uncertainty budget* :math:`b`, you may add the bound :math:`\lVert t \rVert_1 \leq b` (via ``picos.Norm(t, 1) <= b``) which can be represented in conic form as .. math:: &\exists v \in \mathbb{R}^{\operatorname{dim}(t)} : -t \leq v \land t \leq v \land \mathbf{1}^Tv \leq b \\ \Longleftrightarrow~ &\exists v \in \mathbb{R}^{\operatorname{dim}(t)} : \begin{pmatrix} v + t \\ v - t \\ b - \mathbf{1}^Tv \end{pmatrix} \in \mathbb{R}_{\geq 0}^{2\operatorname{dim}(t) + 1}. The auxiliary variable :math:`v` then plays the role of (a slice of) :math:`u` in the formal definition of :math:`\Theta`. When you are done adding bounds, you can obtain :math:`\theta` using :meth:`compile`. :raises RuntimeError: If the set was already compiled. """ self._forbid_compiled() vars = constraint.variables if len(vars) != 1 or next(iter(vars)) != self._element: raise ValueError( "The constraint {} does not bound (only) the perturbation " "element {}.".format(constraint, self._element.string)) self._bounds.add_constraint(constraint)
[docs] def compile(self, validate_feasibility=False): r"""Compile the set and return :math:`\theta`. Internally, this computes the matrices :math:`A` and :math:`B`, the vector :math:`c` and the (product) cone :math:`K`. :param bool validate_feasibility: Whether to solve the feasibility problem associated with the bounds on :math:`t` to verify that :math:`\Theta` is nonempty. :returns: An instance of :class:`~.perturbation.Perturbation`. :raises RuntimeError: If the set was already compiled. :raises TypeError: If the bound constraints could not be put into conic form. :raises ValueError: If :math:`\Theta` could not be verified to be nonempty (needs ``validate_feasibility=True``). """ from ...constraints import DummyConstraint from ...modeling import NoStrategyFound self._forbid_compiled() # If no bounds are given, add a DummyConstraint to mark t free. if not self._bounds.constraints: self._bounds.add_constraint(DummyConstraint(self._element)) # Transform all bounds into conic form. try: conic_bounds = self._bounds.conic_form except NoStrategyFound: raise TypeError("Could not find a conic representation for all of " "the bounds on {}.".format(self._element.string)) # Validate bound feasibility if requested. if validate_feasibility: from ...modeling.solution import SS_OPTIMAL, SS_FEASIBLE solution = conic_bounds.solve( primals=None, apply_solution=False, **settings.INTERNAL_OPTIONS) status = solution.primalStatus if status not in (SS_OPTIMAL, SS_FEASIBLE): raise ValueError("Could not verify that the bounds on {} are " "feasible: The solver {} reports a primal solution state of" " {} for the associated feasibility problem.".format( self._element.string, solution.solver, status)) # Form a virtual variable u from the auxiliary variables. u = [var for var in conic_bounds.variables.values() if var is not self._element] # Convert auxiliary variable bounds to additional affine constraints. conic_bounds.add_list_of_constraints([ var.bound_constraint for var in u if var.bound_constraint]) # Reformulate bounds with respect to a single product cone K. A, B, c, K = [], [], [], [] for constraint in conic_bounds.constraints.values(): member, cone = constraint.conic_membership_form if self._element in member._linear_coefs: A.append(member._linear_coefs[self._element]) else: A.append(cvxopt.spmatrix( [], [], [], size=(cone.dim, self._element.dim))) if u: B.append(cvxopt_hcat([ member._linear_coefs[var] if var in member._linear_coefs else cvxopt.spmatrix([], [], [], size=(cone.dim, var.dim)) for var in u])) c.append(member._constant_coef) K.append(cone) self._A = Constant("A", cvxopt_vcat(A)) self._B = Constant("B", cvxopt_vcat(B)) if u else None self._c = Constant("c", cvxopt_vcat(c)) self._K = ProductCone(*K) # Store the bounds in conic form to speed up worst_case. self._conic_bounds = conic_bounds self._compiled = True return self._parameter
@property def distributional(self): """Implement for :class:`~.perturbation.PerturbationUniverse`.""" return False @property def parameter(self): r"""The perturbation parameter :math:`\theta` living in the set. This is the object returned by :meth:`compile`. :raises RuntimeError: If the set has not been compiled. """ self._require_compiled() return self._parameter @property def A(self): r"""The compiled matrix :math:`A`. :raises RuntimeError: If the set has not been compiled. """ self._require_compiled() return self._A @property def B(self): r"""The compiled matrix :math:`B` or :obj:`None` if :math:`l = 0`. :raises RuntimeError: If the set has not been compiled. """ self._require_compiled() return self._B @property def c(self): r"""The compiled vector :math:`c`. :raises RuntimeError: If the set has not been compiled. """ self._require_compiled() return self._c @property def K(self): r"""The compiled (product) cone :math:`K`. :raises RuntimeError: If the set has not been compiled. """ self._require_compiled() return self._K
[docs] def worst_case(self, scalar, direction): """Implement for :class:`~.perturbation.PerturbationUniverse`.""" from ...modeling import SolutionFailure self._require_compiled() self._check_worst_case_argument_scalar(scalar) self._check_worst_case_argument_direction(direction) p = self._parameter P = self._conic_bounds.copy() x = P.variables[] f = scalar.replace_mutables({p: x}) self._check_worst_case_f_and_x(f, x) if (direction == "min" and not f.convex) \ or (direction == "max" and not f.concave): raise IntractableWorstCase("PICOS refuses to compute {}({}) for {} " "as this is a nonconvex problem.".format(direction, f.string, glyphs.element(, self))) P.set_objective(direction, f if direction != "find" else None) try: P.solve(**settings.INTERNAL_OPTIONS) return f.safe_value, x.safe_value except (SolutionFailure, NotValued) as error: raise RuntimeError("Failed to compute {}({}) for {}: {}".format( direction, f.string, glyphs.element(x.string, self), error))
@cached_property def unit_ball_form(self): """A recipe to repose from ellipsoidal to unit norm ball uncertainty. If the set is :attr:`ellipsoidal`, then this is a pair ``(U, M)`` where ``U`` is a :class:`~.pert_conic.UnitBallPerturbationSet` and ``M`` is a dictionary mapping the old :attr:`parameter` to an affine expression of the new parameter that can represent the old parameter in an expression (see :meth:`~.expression.Expression.replace_mutables`). The mapping ``M`` is empty if and only if the perturbation set is already an instance of :class:`~.pert_conic.UnitBallPerturbationSet`. If the uncertainty set is not ellipsoidal, then this is :obj:`None`. See also :attr:`SOCConstraint.unit_ball_form <.con_soc.SOCConstraint.unit_ball_form>`. :Example: >>> from picos import Problem, RealVariable, sum >>> from picos.uncertain import ConicPerturbationSet >>> # Create a conic perturbation set and a refinement recipe. >>> T = ConicPerturbationSet("t", (2, 2)) >>> T.bound(abs(([[1, 2], [3, 4]] ^ T.element) + 1) <= 10) >>> t = T.compile() >>> U, mapping = T.unit_ball_form >>> print(U) {t' : ‖t'‖ ≤ 1} >>> print(mapping) {<2×2 Perturbation: t>: <2×2 Uncertain Affine Expression: t(t')>} >>> # Define and solve a conically uncertain LP. >>> X = RealVariable("X", (2, 2)) >>> P = Problem() >>> P.set_objective("max", sum(X)) >>> _ = P.add_constraint(X + 2*t <= 10) >>> print(repr(P.parameters["t"].universe)) <2×2 Conic Perturbation Set: {t : ‖[2×2]⊙t + [1]‖ ≤ 10}> >>> _ = P.solve(solver="cvxopt") >>> print(X) [-8.00e+00 1.00e+00] [ 4.00e+00 5.50e+00] >>> # Refine the problem to a unit ball uncertain LP. >>> Q = Problem() >>> Q.set_objective("max", sum(X)) >>> _ = Q.add_constraint(X + 2*mapping[t] <= 10) >>> print(repr(Q.parameters["t'"].universe)) <2×2 Unit Ball Perturbation Set: {t' : ‖t'‖ ≤ 1}> >>> _ = Q.solve(solver="cvxopt") >>> print(X) [-8.00e+00 1.00e+00] [ 4.00e+00 5.50e+00] """ self._require_compiled() if self._B is not None: return None K = self._K.refined if not isinstance(K, SecondOrderCone): return None C = self._A*self._element.vec + self._c << K assert isinstance(C, SOCConstraint) try: X, aff_y, y, _ = C.unit_ball_form except ValueError: return None assert X is self._element U = UnitBallPerturbationSet("{}'".format(, X.shape) u = U.parameter replacement = UncertainAffineExpression("{}({})".format(,, X.shape, {(u,): aff_y._linear_coefs[y], (): aff_y._constant_coef}) return U, {self._parameter: replacement} @property def ellipsoidal(self): """Whether the perturbation set is an ellipsoid. If this is true, then a :attr:`unit_ball_form` is available. """ return bool(self.unit_ball_form)
[docs]class UnitBallPerturbationSet(ConicPerturbationSet): r"""Represents perturbation in an Euclidean or Frobenius unit norm ball. This is a :class:`~.pert_conic.ConicPerturbationSet` with fixed form .. math:: \{t \in \mathbb{R}^{m \times n} \mid \lVert t \rVert_F \leq 1\}. After initialization, you can obtain the parameter using :attr:`~.pert_conic.ConicPerturbationSet.parameter`. """
[docs] def __init__(self, parameter_name, shape=(1, 1)): """See :meth:`ConicPerturbationSet.__init__`.""" ConicPerturbationSet.__init__(self, parameter_name, shape) self.bound(self._element << Ball()) self.compile()
@classmethod def _get_type_string_base(cls): return "Unit Ball Perturbation Set" @property def unit_ball_form(self): """Overwrite :attr:`ConicPerturbationSet.unit_ball_form`.""" return self, {}
# -------------------------------------- __all__ = api_end(_API_START, globals())