Source code for picos.expressions.uncertain.pert_moment

# ------------------------------------------------------------------------------
# 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 MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE.  See the 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:`MomentAmbiguitySet`."""

from collections import namedtuple

from ... import glyphs
from ...apidoc import api_end, api_start
from ..cone_trivial import TheField
from import cvxopt_hpd, cvxopt_inverse, load_shape
from ..exp_affine import AffineExpression
from ..expression import Expression
from ..set_ellipsoid import Ellipsoid
from .perturbation import Perturbation, PerturbationUniverse

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

[docs]class MomentAmbiguitySet(PerturbationUniverse): r"""A moment-uncertain description of a random perturbation parameter. :Model of uncertainty: As a distributional ambiguity set, an instance :math:`\mathcal{P}` of this class 1. represents a safety region for a partially known (ambiguous) probability distribution :math:`\Xi \in \mathcal{P}` and 2. provides a random, ambiguously distributed perturbation parameter :math:`\xi \sim \Xi` that can be used to define worst-case-expectation expressions of the form .. math:: \mathop{(\max\;\textit{or}\;\min)}_{\Xi \in \mathcal{P}} \mathbb{E}_\Xi[f(x, \xi)] for a selection of functions :math:`f` and a decision variable :math:`x`. :Definition: Formally, this class can describe ambiguity sets of the form .. math:: \mathcal{P} = \left\{ \Xi \in \mathcal{M} ~\middle|~ \begin{aligned} \mathbb{P}(\xi \in \mathcal{S}) &= 1, \\ \left( \mathbb{E}[\xi] - \mu \right)^T \Sigma^{-1} \left( \mathbb{E}[\xi] - \mu \right) &\leq \alpha, \\ \mathbb{E}[(\xi - \mu)(\xi - \mu)^T] &\preceq \beta \Sigma \end{aligned} \right\} where 1. :math:`\mathcal{M}` is the set of all Borel probability measures on :math:`\mathbb{R}^n`for some :math:`n \in \mathbb{Z}_{\geq 1}`, 2. the sample space :math:`\mathcal{S} \subseteq \mathbb{R}^n` bounds the support of :math:`\Xi` and may be given as either :math:`\mathbb{R}^n` or as an :math:`n`-dimensional ellipsoid, 3. :math:`\mu \in \mathbb{R}^n` and :math:`\Sigma \in \mathbb{S}^n` with :math:`\Sigma \succ 0` are the **nominal** mean and covariance matrix of :math:`\Xi`, respectively, and 4. :math:`\alpha \geq 0` and :math:`\beta \geq 1` are meta-parameters bounding the uncertainty concerning the mean and the covariance matrix, respectively. Unless :math:`\mathcal{S} = \mathbb{R}^n`, this class can also represent the limit cases of :math:`\alpha \to \infty` and :math:`\beta \to \infty` denoting an unbounded mean and covariance matrix, respectively. .. note:: While :math:`\alpha = 0` denotes that the nominal mean :math:`\mu` is certain, there is a subtle difference between setting :math:`\beta = 1` on the one hand and assuming a certain form for the covariance matrix on the other hand: In the former case, the worst case covariance matrix may be Lowener smaller than the nominal one. Setting a lower bound on the covarianve matrix is computationally difficult and not supported. :Supported functions: 1. A squared norm :math:`f(x, \xi) = \lVert A(x, \xi) \rVert_F^2` where :math:`A` is biaffine in :math:`x` and :math:`\xi`. This can be written as ``abs(A)**2`` in Python. 2. A convex piecewise linear function :math:`f(x, \xi) = max_{i=1}^k a_i(x, \xi)` where :math:`a` is biaffine in :math:`x` and :math:`\xi` for all :math:`i \in [k]`. This can be written as ``picos.max([a_1, ..., a_k])`` in Python. 3. A concave piecewise linear function :math:`f(x, \xi) = min_{i=1}^k a_i(x, \xi)` where :math:`a` is biaffine in :math:`x` and :math:`\xi` for all :math:`i \in [k]`. This can be written as ``picos.min([a_1, ..., a_k])`` in Python. :Example: We show that for unbounded mean and support (i.e. :math:`\alpha \to \infty` and :math:`\beta \to \infty`) and for a finite samples space :math:`S`, this distributional ambiguity set can be used in the context of classical (non-distributional) robust optimization applied to least squares problems. >>> from picos import Constant, diag, Ellipsoid, Problem, RealVariable >>> from picos.uncertain import ConicPerturbationSet, MomentAmbiguitySet >>> import numpy >>> numpy.random.seed(1) >>> # Setup data and variables of the nominal least squares problem. >>> n = 3 >>> A = Constant("A", numpy.random.random((n, n))) >>> b = Constant("b", numpy.random.random(n)) >>> x = RealVariable("x", n) >>> # Setup an ellipsoid S bounding the uncertainty in both models. >>> N = n**2 >>> S = Ellipsoid(N, diag(range(1, N + 1)), range(N)) >>> # Define and solve both robust models. >>> U1 = ConicPerturbationSet.from_constraints( ... "Y", RealVariable("Y", N) << S) >>> U2 = MomentAmbiguitySet("Z", N, alpha=None, beta=None, sample_space=S) >>> Y = U1.parameter.reshaped((n, n)) >>> Z = U2.parameter.reshaped((n, n)) >>> P1, P2 = Problem(), Problem() >>> P1.objective = "min", abs((A + Y)*x - b) >>> P2.objective = "min", abs((A + Z)*x - b)**2 >>> _ = P1.solve(solver="cvxopt") >>> x1 = ~x # Save current value of x as a constant PICOS expression x1. >>> _ = P2.solve(solver="cvxopt") >>> x2 = ~x >>> # Verify that both models yield the same robust regression vector. >>> x1.equals(x2, relTol=1e-4) True """
[docs] def __init__(self, parameter_name, shape, nominal_mean=0, nominal_covariance="I", alpha=0, beta=1, sample_space=None): r"""Create a :class:`MomentAmbiguitySet`. :param str parameter_name: Name of the random parameter :math:`\xi`. :param shape: Shape of :math:`\xi`. Must characterize a column vector. If :obj:`None`, then the shape is inferred from the nominal mean. :type shape: anything recognized by :func:`` :param nominal_mean: The nominal mean :math:`\mu` of the ambiguous distribution :math:`\Xi`. :type nominal_mean: anything recognized by :func:`` :param nominal_covariance: The nominal covariance matrix :math:`\Sigma` of the ambiguous distribution :math:`\Xi`. :type nominal_covariance: anything recognized by :func:`` :param float alpha: The parameter :math:`\alpha \geq 0` bounding the uncertain mean. The values :obj:`None` and ``float("inf")`` denote an unbounded mean. :param float beta: The parameter :math:`\beta \geq 1` bounding the uncertain covariance matrix. The values :obj:`None` and ``float("inf")`` denote unbounded covariances. :param sample_space: The sample space :math:`\mathcal{S}`. If this is :obj:`None` or an instance of :class:`~picos.TheField` (i.e. :math:`\mathbb{R}^n`), then the support of :math:`\Xi` is unbounded. If this is an :math:`n`-dimensional instance of :class:`~picos.Ellipsoid`, then the support of :math:`\Xi` is a subset of that ellipsoid. :type sample_space: :obj:`None` or :class:`~picos.TheField` or :class:`~picos.Ellipsoid` """ # Load the dimensionality. shape = load_shape(shape) if shape[1] != 1: raise ValueError( "The perturbation parameter must be a column vector; a shape of" " {} is not supported.".format(glyphs.shape(shape))) self._dim = n = shape[0] # Load the nominal mean. if not isinstance(nominal_mean, Expression): nominal_mean = AffineExpression.from_constant( nominal_mean, shape=shape) else: nominal_mean = nominal_mean.refined if not isinstance(nominal_mean, AffineExpression) \ or not nominal_mean.constant: raise TypeError("The nominal mean must be a real constant.") if nominal_mean.shape != shape: raise TypeError("The nominal mean must be a {}-dimensional " "column vector, got {} instead." .format(n, glyphs.shape(nominal_mean.shape))) self._mean = nominal_mean # Load the nominal covariance matrix. if not isinstance(nominal_covariance, Expression): nominal_covariance = AffineExpression.from_constant( nominal_covariance, shape=(n, n)) else: nominal_covariance = nominal_covariance.refined if not isinstance(nominal_covariance, AffineExpression) \ or not nominal_covariance.constant: raise TypeError( "The nominal covarance matrix must be a real constant.") if nominal_covariance.shape != (n, n): raise TypeError("The nominal covariance matrix must be a {} " "matrix, got {} instead.".format(glyphs.shape((n, n)), glyphs.shape(nominal_covariance.shape))) cov_value = nominal_covariance.safe_value_as_matrix if not cvxopt_hpd(cov_value): raise ValueError("The nominal covariance matrix is not symmetric " "positive definite.") self._cov = nominal_covariance # Compute the inverse of the nominal covariance matrix. self._cov_inv = AffineExpression.from_constant( cvxopt_inverse(cov_value), shape=(n, n)) # Load alpha. if alpha in (None, float("inf")): self._alpha = None else: self._alpha = float(alpha) if self._alpha < 0: raise ValueError("The parameter alpha must be nonnegative.") # Load beta. if beta in (None, float("inf")): self._beta = None else: self._beta = float(beta) if self._beta < 1: raise ValueError("The parameter beta must be at least one.") # Load the sample space. if sample_space is None: self._ss = None elif isinstance(sample_space, TheField): if sample_space.dim is not None and sample_space.dim != self._dim: raise TypeError("Invalid sample space dimension.") self._ss = None elif isinstance(sample_space, Ellipsoid): if sample_space.dim != self._dim: raise TypeError("Invalid sample space dimension.") self._ss = sample_space else: raise TypeError("Invalid sample space type.") # Limit the permissible combinations of unboundedness. if self._ss is None and (self._alpha is None or self._beta is None): raise ValueError("If the support is unbounded, both the mean and " "the covariance matrix must be bounded.") # Create the perturbation parameter. self._parameter = Perturbation(self, parameter_name, shape)
@property def dim(self): """The dimension :math:`n` of the sample space.""" return self._dim @property def nominal_mean(self): r"""The nominal mean :math:`\mu` of the ambiguous distribution.""" return self._mean @property def nominal_covariance(self): r"""The nominal covariance matrix :math:`\Sigma`.""" return self._cov @property def alpha(self): r"""The parameter :math:`\alpha`. A value of :obj:`None` denotes :math:`\alpha \to \infty`. """ return self._alpha @property def beta(self): r"""The parameter :math:`\beta`. A value of :obj:`None` denotes :math:`\beta \to \infty`. """ return self._beta @property def sample_space(self): r"""The sample space (bound on support) :math:`\mathcal{S}`. A value of :obj:`None` means :math:`\mathcal{S} = \mathbb{R}^n`. """ return self._ss Subtype = namedtuple("Subtype", ( "dim", "bounded_mean", "bounded_covariance", "bounded_support")) def _subtype(self): return self.Subtype( self._dim, self._alpha is not None, self._beta is not None, self._ss is not None) def __str__(self): return "MAS(mu={}, cov={}, a={}, b={}, S={})".format( self._mean.string, self._cov.string, self._alpha, self._beta, self._ss.string if self._ss is not None else None) @classmethod def _get_type_string_base(cls): return "Moment Ambiguity Set" def __repr__(self): return glyphs.repr2("{} {}".format(glyphs.shape(self._parameter.shape), self._get_type_string_base()), self.__str__()) @property def distributional(self): """Implement for :class:`~.perturbation.PerturbationUniverse`.""" return True @property def parameter(self): r"""The random perturbation parameter :math:`\xi`.""" return self._parameter
# -------------------------------------- __all__ = api_end(_API_START, globals())