Coverage for picos/expressions/uncertain/pert_moment.py: 75.00%
96 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +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:`MomentAmbiguitySet`."""
21from collections import namedtuple
23from ... import glyphs
24from ...apidoc import api_end, api_start
25from ..cone_trivial import TheField
26from ..data import cvxopt_hpd, cvxopt_inverse, load_shape
27from ..exp_affine import AffineExpression
28from ..expression import Expression
29from ..set_ellipsoid import Ellipsoid
30from .perturbation import Perturbation, PerturbationUniverse
32_API_START = api_start(globals())
33# -------------------------------
36class MomentAmbiguitySet(PerturbationUniverse):
37 r"""A moment-uncertain description of a random perturbation parameter.
39 :Model of uncertainty:
41 As a distributional ambiguity set, an instance :math:`\mathcal{P}` of this
42 class
44 1. represents a safety region for a partially known (ambiguous) probability
45 distribution :math:`\Xi \in \mathcal{P}` and
46 2. provides a random, ambiguously distributed perturbation parameter
47 :math:`\xi \sim \Xi` that can be used to define worst-case-expectation
48 expressions of the form
50 .. math::
52 \mathop{(\max\;\textit{or}\;\min)}_{\Xi \in \mathcal{P}}
53 \mathbb{E}_\Xi[f(x, \xi)]
55 for a selection of functions :math:`f` and a decision variable :math:`x`.
57 :Definition:
59 Formally, this class can describe ambiguity sets of the form
61 .. math::
63 \mathcal{P} = \left\{
64 \Xi \in \mathcal{M} ~\middle|~
65 \begin{aligned}
66 \mathbb{P}(\xi \in \mathcal{S}) &= 1, \\
67 \left( \mathbb{E}[\xi] - \mu \right)^T \Sigma^{-1}
68 \left( \mathbb{E}[\xi] - \mu \right) &\leq \alpha, \\
69 \mathbb{E}[(\xi - \mu)(\xi - \mu)^T] &\preceq \beta \Sigma
70 \end{aligned}
71 \right\}
73 where
75 1. :math:`\mathcal{M}` is the set of all Borel probability measures on
76 :math:`\mathbb{R}^n`for some :math:`n \in \mathbb{Z}_{\geq 1}`,
77 2. the sample space :math:`\mathcal{S} \subseteq \mathbb{R}^n` bounds the
78 support of :math:`\Xi` and may be given as either :math:`\mathbb{R}^n`
79 or as an :math:`n`-dimensional ellipsoid,
80 3. :math:`\mu \in \mathbb{R}^n` and
81 :math:`\Sigma \in \mathbb{S}^n` with :math:`\Sigma \succ 0` are the
82 **nominal** mean and covariance matrix of :math:`\Xi`, respectively, and
83 4. :math:`\alpha \geq 0` and :math:`\beta \geq 1` are meta-parameters
84 bounding the uncertainty concerning the mean and the covariance matrix,
85 respectively.
87 Unless :math:`\mathcal{S} = \mathbb{R}^n`, this class can also represent the
88 limit cases of :math:`\alpha \to \infty` and :math:`\beta \to \infty`
89 denoting an unbounded mean and covariance matrix, respectively.
91 .. note::
93 While :math:`\alpha = 0` denotes that the nominal mean :math:`\mu` is
94 certain, there is a subtle difference between setting :math:`\beta = 1`
95 on the one hand and assuming a certain form for the covariance matrix on
96 the other hand: In the former case, the worst case covariance matrix may
97 be Lowener smaller than the nominal one. Setting a lower bound on the
98 covarianve matrix is computationally difficult and not supported.
100 :Supported functions:
102 1. A squared norm :math:`f(x, \xi) = \lVert A(x, \xi) \rVert_F^2` where
103 :math:`A` is biaffine in :math:`x` and :math:`\xi`. This can be written
104 as ``abs(A)**2`` in Python.
105 2. A convex piecewise linear function :math:`f(x, \xi) = max_{i=1}^k a_i(x,
106 \xi)` where :math:`a` is biaffine in :math:`x` and :math:`\xi` for all
107 :math:`i \in [k]`. This can be written as ``picos.max([a_1, ..., a_k])``
108 in Python.
109 3. A concave piecewise linear function :math:`f(x, \xi) = min_{i=1}^k a_i(x,
110 \xi)` where :math:`a` is biaffine in :math:`x` and :math:`\xi` for all
111 :math:`i \in [k]`. This can be written as ``picos.min([a_1, ..., a_k])``
112 in Python.
114 :Example:
116 We show that for unbounded mean and support (i.e. :math:`\alpha \to \infty`
117 and :math:`\beta \to \infty`) and for a finite samples space :math:`S`, this
118 distributional ambiguity set can be used in the context of classical
119 (non-distributional) robust optimization applied to least squares problems.
121 >>> from picos import Constant, diag, Ellipsoid, Problem, RealVariable
122 >>> from picos.uncertain import ConicPerturbationSet, MomentAmbiguitySet
123 >>> import numpy
124 >>> numpy.random.seed(1)
125 >>> # Setup data and variables of the nominal least squares problem.
126 >>> n = 3
127 >>> A = Constant("A", numpy.random.random((n, n)))
128 >>> b = Constant("b", numpy.random.random(n))
129 >>> x = RealVariable("x", n)
130 >>> # Setup an ellipsoid S bounding the uncertainty in both models.
131 >>> N = n**2
132 >>> S = Ellipsoid(N, diag(range(1, N + 1)), range(N))
133 >>> # Define and solve both robust models.
134 >>> U1 = ConicPerturbationSet.from_constraints(
135 ... "Y", RealVariable("Y", N) << S)
136 >>> U2 = MomentAmbiguitySet("Z", N, alpha=None, beta=None, sample_space=S)
137 >>> Y = U1.parameter.reshaped((n, n))
138 >>> Z = U2.parameter.reshaped((n, n))
139 >>> P1, P2 = Problem(), Problem()
140 >>> P1.objective = "min", abs((A + Y)*x - b)
141 >>> P2.objective = "min", abs((A + Z)*x - b)**2
142 >>> _ = P1.solve(solver="cvxopt")
143 >>> x1 = ~x # Save current value of x as a constant PICOS expression x1.
144 >>> _ = P2.solve(solver="cvxopt")
145 >>> x2 = ~x
146 >>> # Verify that both models yield the same robust regression vector.
147 >>> x1.equals(x2, relTol=1e-4)
148 True
149 """
151 def __init__(self, parameter_name, shape, nominal_mean=0,
152 nominal_covariance="I", alpha=0, beta=1, sample_space=None):
153 r"""Create a :class:`MomentAmbiguitySet`.
155 :param str parameter_name:
156 Name of the random parameter :math:`\xi`.
158 :param shape:
159 Shape of :math:`\xi`. Must characterize a column vector. If
160 :obj:`None`, then the shape is inferred from the nominal mean.
161 :type shape:
162 anything recognized by :func:`~picos.expressions.data.load_shape`
164 :param nominal_mean:
165 The nominal mean :math:`\mu` of the ambiguous distribution
166 :math:`\Xi`.
167 :type nominal_mean:
168 anything recognized by :func:`~picos.expressions.data.load_data`
170 :param nominal_covariance:
171 The nominal covariance matrix :math:`\Sigma` of the ambiguous
172 distribution :math:`\Xi`.
173 :type nominal_covariance:
174 anything recognized by :func:`~picos.expressions.data.load_data`
176 :param float alpha:
177 The parameter :math:`\alpha \geq 0` bounding the uncertain mean.
178 The values :obj:`None` and ``float("inf")`` denote an unbounded
179 mean.
181 :param float beta:
182 The parameter :math:`\beta \geq 1` bounding the uncertain covariance
183 matrix. The values :obj:`None` and ``float("inf")`` denote unbounded
184 covariances.
186 :param sample_space:
187 The sample space :math:`\mathcal{S}`. If this is :obj:`None` or an
188 instance of :class:`~picos.TheField` (i.e. :math:`\mathbb{R}^n`),
189 then the support of :math:`\Xi` is unbounded. If this is an
190 :math:`n`-dimensional instance of :class:`~picos.Ellipsoid`, then
191 the support of :math:`\Xi` is a subset of that ellipsoid.
192 :type sample_space:
193 :obj:`None` or :class:`~picos.TheField` or :class:`~picos.Ellipsoid`
194 """
195 # Load the dimensionality.
196 shape = load_shape(shape)
197 if shape[1] != 1:
198 raise ValueError(
199 "The perturbation parameter must be a column vector; a shape of"
200 " {} is not supported.".format(glyphs.shape(shape)))
201 self._dim = n = shape[0]
203 # Load the nominal mean.
204 if not isinstance(nominal_mean, Expression):
205 nominal_mean = AffineExpression.from_constant(
206 nominal_mean, shape=shape)
207 else:
208 nominal_mean = nominal_mean.refined
210 if not isinstance(nominal_mean, AffineExpression) \
211 or not nominal_mean.constant:
212 raise TypeError("The nominal mean must be a real constant.")
214 if nominal_mean.shape != shape:
215 raise TypeError("The nominal mean must be a {}-dimensional "
216 "column vector, got {} instead."
217 .format(n, glyphs.shape(nominal_mean.shape)))
219 self._mean = nominal_mean
221 # Load the nominal covariance matrix.
222 if not isinstance(nominal_covariance, Expression):
223 nominal_covariance = AffineExpression.from_constant(
224 nominal_covariance, shape=(n, n))
225 else:
226 nominal_covariance = nominal_covariance.refined
228 if not isinstance(nominal_covariance, AffineExpression) \
229 or not nominal_covariance.constant:
230 raise TypeError(
231 "The nominal covarance matrix must be a real constant.")
233 if nominal_covariance.shape != (n, n):
234 raise TypeError("The nominal covariance matrix must be a {} "
235 "matrix, got {} instead.".format(glyphs.shape((n, n)),
236 glyphs.shape(nominal_covariance.shape)))
238 cov_value = nominal_covariance.safe_value_as_matrix
240 if not cvxopt_hpd(cov_value):
241 raise ValueError("The nominal covariance matrix is not symmetric "
242 "positive definite.")
244 self._cov = nominal_covariance
246 # Compute the inverse of the nominal covariance matrix.
247 self._cov_inv = AffineExpression.from_constant(
248 cvxopt_inverse(cov_value), shape=(n, n))
250 # Load alpha.
251 if alpha in (None, float("inf")):
252 self._alpha = None
253 else:
254 self._alpha = float(alpha)
255 if self._alpha < 0:
256 raise ValueError("The parameter alpha must be nonnegative.")
258 # Load beta.
259 if beta in (None, float("inf")):
260 self._beta = None
261 else:
262 self._beta = float(beta)
263 if self._beta < 1:
264 raise ValueError("The parameter beta must be at least one.")
266 # Load the sample space.
267 if sample_space is None:
268 self._ss = None
269 elif isinstance(sample_space, TheField):
270 if sample_space.dim is not None and sample_space.dim != self._dim:
271 raise TypeError("Invalid sample space dimension.")
273 self._ss = None
274 elif isinstance(sample_space, Ellipsoid):
275 if sample_space.dim != self._dim:
276 raise TypeError("Invalid sample space dimension.")
278 self._ss = sample_space
279 else:
280 raise TypeError("Invalid sample space type.")
282 # Limit the permissible combinations of unboundedness.
283 if self._ss is None and (self._alpha is None or self._beta is None):
284 raise ValueError("If the support is unbounded, both the mean and "
285 "the covariance matrix must be bounded.")
287 # Create the perturbation parameter.
288 self._parameter = Perturbation(self, parameter_name, shape)
290 @property
291 def dim(self):
292 """The dimension :math:`n` of the sample space."""
293 return self._dim
295 @property
296 def nominal_mean(self):
297 r"""The nominal mean :math:`\mu` of the ambiguous distribution."""
298 return self._mean
300 @property
301 def nominal_covariance(self):
302 r"""The nominal covariance matrix :math:`\Sigma`."""
303 return self._cov
305 @property
306 def alpha(self):
307 r"""The parameter :math:`\alpha`.
309 A value of :obj:`None` denotes :math:`\alpha \to \infty`.
310 """
311 return self._alpha
313 @property
314 def beta(self):
315 r"""The parameter :math:`\beta`.
317 A value of :obj:`None` denotes :math:`\beta \to \infty`.
318 """
319 return self._beta
321 @property
322 def sample_space(self):
323 r"""The sample space (bound on support) :math:`\mathcal{S}`.
325 A value of :obj:`None` means :math:`\mathcal{S} = \mathbb{R}^n`.
326 """
327 return self._ss
329 Subtype = namedtuple("Subtype", (
330 "dim",
331 "bounded_mean",
332 "bounded_covariance",
333 "bounded_support"))
335 def _subtype(self):
336 return self.Subtype(
337 self._dim,
338 self._alpha is not None,
339 self._beta is not None,
340 self._ss is not None)
342 def __str__(self):
343 return "MAS(mu={}, cov={}, a={}, b={}, S={})".format(
344 self._mean.string,
345 self._cov.string,
346 self._alpha,
347 self._beta,
348 self._ss.string if self._ss is not None else None)
350 @classmethod
351 def _get_type_string_base(cls):
352 return "Moment Ambiguity Set"
354 def __repr__(self):
355 return glyphs.repr2("{} {}".format(glyphs.shape(self._parameter.shape),
356 self._get_type_string_base()), self.__str__())
358 @property
359 def distributional(self):
360 """Implement for :class:`~.perturbation.PerturbationUniverse`."""
361 return True
363 @property
364 def parameter(self):
365 r"""The random perturbation parameter :math:`\xi`."""
366 return self._parameter
369# --------------------------------------
370__all__ = api_end(_API_START, globals())