picos.uncertain

Tools for handling uncertain data.

Outline

Classes

ConicPerturbationSet

A conic description of a Perturbation.

MomentAmbiguitySet

A moment-uncertain description of a random perturbation parameter.

ScenarioPerturbationSet

A scenario description of a Perturbation.

UnitBallPerturbationSet

Represents perturbation in an Euclidean or Frobenius unit norm ball.

WassersteinAmbiguitySet

A wasserstein ambiguity set centered at a discrete distribution.

Exceptions

IntractableWorstCase

Computing a worst-case (expected) value is hard and not supported.

Classes

ConicPerturbationSet

class picos.uncertain.ConicPerturbationSet(parameter_name, shape=(1, 1))[source]

Bases: picos.expressions.uncertain.perturbation.PerturbationUniverse

A conic description of a Perturbation.

Definition

An instance \Theta of this class defines a perturbation parameter

\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 m, n \in \mathbb{Z}_{\geq 1}, l \in \mathbb{Z}_{\geq 0}, A \in \mathbb{R}^{k \times mn}, B \in \mathbb{R}^{k \times l}, c \in \mathbb{R}^k and K \subseteq \mathbb{R}^k is a (product) cone for some k \in \mathbb{Z}_{\geq 1}.

Usage

Obtaining \theta is done in a number of steps:

  1. Create an instance of this class (see __init__).

  2. Access element to obtain a regular, fresh RealVariable representing t.

  3. Define \Theta through any number of regular PICOS constraints that depend only on t and that have a conic representation by passing the constraints to bound.

  4. Call compile to obtain a handle to the Perturbation \theta.

  5. You can now use \theta to build instances of UncertainAffineExpression and derived constraint types.

It is best practice to assign t to a Python variable and overwrite that variable with \theta on compilation.

Alternatively, you can obtain a compiled \Theta from the factory method from_constraints and access \theta via 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>
__init__(parameter_name, shape=(1, 1))[source]

Create a ConicPerturbationSet.

Parameters
  • parameter_name (str) – Name of the parameter that lives in the set.

  • shape (int or tuple or list) – The shape of a vector or matrix perturbation.

bound(constraint)[source]

Add a constraint that bounds 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 b, you may add the bound \lVert t \rVert_1 \leq b (via picos.Norm(t, 1) <= b) which can be represented in conic form as

&\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 v then plays the role of (a slice of) u in the formal definition of \Theta.

When you are done adding bounds, you can obtain \theta using compile.

Raises

RuntimeError – If the set was already compiled.

compile(validate_feasibility=False)[source]

Compile the set and return \theta.

Internally, this computes the matrices A and B, the vector c and the (product) cone K.

Parameters

validate_feasibility (bool) – Whether to solve the feasibility problem associated with the bounds on t to verify that \Theta is nonempty.

Returns

An instance of Perturbation.

Raises
  • RuntimeError – If the set was already compiled.

  • TypeError – If the bound constraints could not be put into conic form.

  • ValueError – If \Theta could not be verified to be nonempty (needs validate_feasibility=True).

classmethod from_constraints(parameter_name, *constraints)[source]

Create a ConicPerturbationSet from constraints.

The constraints must concern a single regular decision variable that plays the role of the element t. This variable is not stored or modified and can be reused in a different context.

Parameters
  • parameter_name (str) – Name of the parameter that lives in the set.

  • 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.

  • 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>
worst_case(scalar, direction)[source]

Implement for PerturbationUniverse.

property A

The compiled matrix A.

Raises

RuntimeError – If the set has not been compiled.

property B

The compiled matrix B or None if l = 0.

Raises

RuntimeError – If the set has not been compiled.

property K

The compiled (product) cone K.

Raises

RuntimeError – If the set has not been compiled.

property c

The compiled vector c.

Raises

RuntimeError – If the set has not been compiled.

property distributional

Implement for PerturbationUniverse.

property element

The perturbation element t describing the set.

This is a regular RealVariable that you can use to create constraints to pass to bound. You can then obtain the “actual” perturbation parameter \theta to use in expressions alongside your decision variaiables using compile.

Warning

If you use this object instead of 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.

property ellipsoidal

Whether the perturbation set is an ellipsoid.

If this is true, then a unit_ball_form is available.

property parameter

The perturbation parameter \theta living in the set.

This is the object returned by compile.

Raises

RuntimeError – If the set has not been compiled.

property unit_ball_form

A recipe to repose from ellipsoidal to unit norm ball uncertainty.

If the set is ellipsoidal, then this is a pair (U, M) where U is a UnitBallPerturbationSet and M is a dictionary mapping the old parameter to an affine expression of the new parameter that can represent the old parameter in an expression (see replace_mutables). The mapping M is empty if and only if the perturbation set is already an instance of UnitBallPerturbationSet.

If the uncertainty set is not ellipsoidal, then this is None.

See also 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]

MomentAmbiguitySet

class picos.uncertain.MomentAmbiguitySet(parameter_name, shape, nominal_mean=0, nominal_covariance='I', alpha=0, beta=1, sample_space=None)[source]

Bases: picos.expressions.uncertain.perturbation.PerturbationUniverse

A moment-uncertain description of a random perturbation parameter.

Model of uncertainty

As a distributional ambiguity set, an instance \mathcal{P} of this class

  1. represents a safety region for a partially known (ambiguous) probability distribution \Xi \in \mathcal{P} and

  2. provides a random, ambiguously distributed perturbation parameter \xi \sim \Xi that can be used to define worst-case-expectation expressions of the form

    \mathop{(\max\;\textit{or}\;\min)}_{\Xi \in \mathcal{P}}
\mathbb{E}_\Xi[f(x, \xi)]

    for a selection of functions f and a decision variable x.

Definition

Formally, this class can describe ambiguity sets of the form

\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. \mathcal{M} is the set of all Borel probability measures on \mathbb{R}^n,

  2. the sample space \mathcal{S} \subseteq \mathbb{R}^n bounds the support of \Xi and may be given as either \mathbb{R}^n or as an n-dimensional ellipsoid,

  3. \mu \in \mathbb{R}^n and \Sigma \in \mathbb{S}^n with \Sigma \succ 0 are the nominal mean and covariance matrix of \Xi, respectively, and

  4. \alpha \geq 0 and \beta \geq 1 are meta-parameters bounding the uncertainty concerning the mean and the covariance matrix, respectively.

Unless \mathcal{S} = \mathbb{R}^n, this class can also represent the limit cases of \alpha \to \infty and \beta \to \infty denoting an unbounded mean and covariance matrix, respectively.

Note

While \alpha = 0 denotes that the nominal mean \mu is certain, there is a subtle difference between setting \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 f(x, \xi) = \lVert A(x, \xi) \rVert_F^2 where A is biaffine in x and \xi. This can be written as abs(A)**2 in Python.

  2. A convex piecewise linear function f(x, \xi) = max_{i=1}^k a_i(x,
\xi) where a is biaffine in x and \xi for all i \in [k]. This can be written as picos.max([a_1, ..., a_k]) in Python.

  3. A concave piecewise linear function f(x, \xi) = min_{i=1}^k a_i(x,
\xi) where a is biaffine in x and \xi for all 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. \alpha \to \infty and \beta \to \infty) and for a finite samples space 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
__init__(parameter_name, shape, nominal_mean=0, nominal_covariance='I', alpha=0, beta=1, sample_space=None)[source]

Create a MomentAmbiguitySet.

Parameters
  • parameter_name (str) – Name of the random parameter \xi.

  • shape (anything recognized by load_shape) – Shape of \xi. Must characterize a column vector. If None, then the shape is inferred from the nominal mean.

  • nominal_mean (anything recognized by load_data) – The nominal mean \mu of the ambiguous distribution \Xi.

  • nominal_covariance (anything recognized by load_data) – The nominal covariance matrix \Sigma of the ambiguous distribution \Xi.

  • alpha (float) – The parameter \alpha \geq 0 bounding the uncertain mean. The values None and float("inf") denote an unbounded mean.

  • beta (float) – The parameter \beta \geq 1 bounding the uncertain covariance matrix. The values None and float("inf") denote unbounded covariances.

  • sample_space (None or TheField or Ellipsoid) – The sample space \mathcal{S}. If this is None or an instance of TheField (i.e. \mathbb{R}^n), then the support of \Xi is unbounded. If this is an n-dimensional instance of Ellipsoid, then the support of \Xi is a subset of that ellipsoid.

property alpha

The parameter \alpha.

A value of None denotes \alpha \to \infty.

property beta

The parameter \beta.

A value of None denotes \beta \to \infty.

property dim

The dimension n of the sample space.

property distributional

Implement for PerturbationUniverse.

property nominal_covariance

The nominal covariance matrix \Sigma.

property nominal_mean

The nominal mean \mu of the ambiguous distribution.

property parameter

The random perturbation parameter \xi.

property sample_space

The sample space (bound on support) \mathcal{S}.

A value of None means \mathcal{S} = \mathbb{R}^n.

ScenarioPerturbationSet

class picos.uncertain.ScenarioPerturbationSet(parameter_name, scenarios, compute_hull=None)[source]

Bases: picos.expressions.uncertain.perturbation.PerturbationUniverse

A scenario description of a Perturbation.

Definition

An instance \Theta of this class defines a perturbation parameter

\theta \in \Theta = \operatorname{conv}(S)

where S \subset \mathbb{R}^{m \times n} is a finite set of scenarios and \operatorname{conv} denotes the convex hull.

Usually, the scenarios are observed or projected realizations of the uncertain data and \theta is used to represent the data directly.

Example

>>> from picos.uncertain import ScenarioPerturbationSet
>>> scenarios = [[1, -1], [1, 1], [-1, -1], [-1, 1], [0, 0]]
>>> S = ScenarioPerturbationSet("s", scenarios, False); S
<2×1 Scenario Perturbation Set: conv({5 2×1 scenarios})>
>>> s = S.parameter; s
<2×1 Perturbation: s>
>>> # Compute largest sum of entries over all points in S.
>>> value, realization = S.worst_case(s[0] + s[1], "max")
>>> round(value,  4)
2.0
>>> print(realization)
[ 1.00e+00]
[ 1.00e+00]
__init__(parameter_name, scenarios, compute_hull=None)[source]

Create a ScenarioPerturbationSet.

Parameters
  • parameter_name (str) – Name of the parameter that lives in the set.

  • scenarios (anything recognized by picos.Samples) – A collection of data points of same shape representing S.

  • compute_hull (bool) – Whether to use SciPy to compute the convex hull of the data points and discard points in the interior. This can speed up the solution process significantly, in particular when the scenarios come from observations and when the data is low-dimensional. On the other hand, when the given scenarios are known to be on the boundary of their convex hull, then disabling this speeds up initialization of the perturbation set. The default value of None means True when SciPy is available and False otherwise.

worst_case(scalar, direction)[source]

Implement for PerturbationUniverse.

property distributional

Implement for PerturbationUniverse.

property parameter

Implement for PerturbationUniverse.

property scenarios

The registered scenarios as a Samples object.

UnitBallPerturbationSet

class picos.uncertain.UnitBallPerturbationSet(parameter_name, shape=(1, 1))[source]

Bases: picos.expressions.uncertain.pert_conic.ConicPerturbationSet

Represents perturbation in an Euclidean or Frobenius unit norm ball.

This is a ConicPerturbationSet with fixed form

\{t \in \mathbb{R}^{m \times n} \mid \lVert t \rVert_F \leq 1\}.

After initialization, you can obtain the parameter using parameter.

__init__(parameter_name, shape=(1, 1))[source]

See ConicPerturbationSet.__init__.

property unit_ball_form

Overwrite ConicPerturbationSet.unit_ball_form.

WassersteinAmbiguitySet

class picos.uncertain.WassersteinAmbiguitySet(parameter_name, p, eps, samples, weights=1)[source]

Bases: picos.expressions.uncertain.perturbation.PerturbationUniverse

A wasserstein ambiguity set centered at a discrete distribution.

Model of uncertainty

As a distributional ambiguity set, an instance \mathcal{P} of this class

  1. represents a safety region for a partially known (ambiguous) probability distribution \Xi \in \mathcal{P} and

  2. provides a random, ambiguously distributed perturbation parameter \xi \sim \Xi that can be used to define worst-case-expectation expressions of the form

    \mathop{(\max\;\textit{or}\;\min)}_{\Xi \in \mathcal{P}}
\mathbb{E}_\Xi[f(x, \xi)]

    for a selection of functions f and a decision variable x.

Definition

Formally, this class can describe discrepancy-based ambiguity sets of the form

\mathcal{P} = \left\{
    \Xi \in \mathcal{M} ~\middle|~
    \operatorname{W}_p(\Xi, \Xi_\text{N}) \leq \epsilon
\right\}

where discrepancy from the discrete nominal distribution

\Xi_\text{N} = \sum_{i = 1}^N w_i \delta_{\xi_{(i)}} \in \mathcal{M}

is measured with respect to the Wasserstein distance of order p \geq 1,

\operatorname{W}_p(\Xi, \Xi') =
    {\left(
        \inf_{\Phi \in \Pi(\Xi, \Xi')}
        \int_{\mathbb{R}^m \times \mathbb{R}^m}
        \lVert \phi - \phi' \rVert^p \;
        \Phi(
            \mathop{}\!\mathrm{d} \phi
            \times
            \mathop{}\!\mathrm{d} \phi')
    \right)}^{\frac{1}{p}},

where

  1. \mathcal{M} is the set of all Borel probability measures on \mathbb{R}^n for some n \in \mathbb{Z}_{\geq 1},

  2. \Pi(\Xi, \Xi') denotes the set of all couplings of \Xi and \Xi',

  3. \xi_{(i)} \in \mathbb{R}^n for all i \in [N] are the N \in \mathbb{Z}_{\geq 1} samples comprising the support of \Xi_\text{N},

  4. w_i \in \mathbb{R}_{\geq 0} are weights denoting the nominal probabilitiy mass at \xi_{(i)} for all i \in [N],

  5. \delta_{\xi_{(i)}} denotes the Dirac delta function with unit mass at \xi_{(i)} for all i \in [N] and where

  6. \epsilon \in \mathbb{R}_{\geq 0} controls the radius of the ambiguity set.

Supported functions

For p = 1:

  1. A convex piecewise linear function f(x, \xi) = max_{i=1}^k a_i(x,
\xi) where a is biaffine in x and \xi for all i \in [k]. This can be written as picos.max([a_1, ..., a_k]) in Python.

  2. A concave piecewise linear function f(x, \xi) = min_{i=1}^k a_i(x,
\xi) where a is biaffine in x and \xi for all i \in [k]. This can be written as picos.min([a_1, ..., a_k]) in Python.

For p = 2:

  1. A squared norm f(x, \xi) = \lVert A(x, \xi) \rVert_F^2 where A is biaffine in x and \xi. This can be written as abs(A)**2 in Python.

__init__(parameter_name, p, eps, samples, weights=1)[source]

Create a WassersteinAmbiguitySet.

Parameters
  • parameter_name (str) – Name of the random parameter \xi.

  • p (float) – The Wasserstein type/order parameter p.

  • eps (float) – The Wasserstein ball radius \epsilon.

  • samples (aynthing recognized by Samples) – The support of the discrete distribution \Xi_\text{D} given as the samples \xi_{(i)}. The original shape of the samples determines the shape of \xi.

  • weights – A vector denoting the nonnegative weight (e.g. frequency or probability) of each sample. Its length must match the number of samples provided. The argument will be normalized such that its entries sum to one. Entries of zero will be dropped alongside their associated sample. The default value of 1 denotes the empirical distribution on the samples.

Warning

Duplicate samples are not detected and can impact performance. If duplicate samples are likely, make sure to detect them and encode their frequency in the weight vector.

property distributional

Implement for PerturbationUniverse.

property eps

The Wasserstein ball radius \epsilon.

property p

The Wasserstein order p.

property parameter

The random perturbation parameter \xi.

property samples

The registered samples as a Samples object.

property weights

The sample weights a constant PICOS vector.

Exceptions

IntractableWorstCase

exception picos.uncertain.IntractableWorstCase[source]

Bases: RuntimeError

Computing a worst-case (expected) value is hard and not supported.

Raised by worst_case and methods that depend on it.

__init__(*args, **kwargs)

Initialize self. See help(type(self)) for accurate signature.

__new__(**kwargs)

Create and return a new object. See help(type) for accurate signature.