Coverage for picos/expressions/uncertain/pert_scenario.py: 88.89%
90 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:`ScenarioPerturbationSet`."""
21from collections import namedtuple
23import cvxopt
25from ... import glyphs, settings
26from ...apidoc import api_end, api_start
27from ..data import cvx2np
28from ..expression import NotValued
29from ..samples import Samples
30from ..variables import RealVariable
31from .perturbation import Perturbation, PerturbationUniverse
32from .uexpression import IntractableWorstCase
34_API_START = api_start(globals())
35# -------------------------------
38class ScenarioPerturbationSet(PerturbationUniverse):
39 r"""A scenario description of a :class:`~.perturbation.Perturbation`.
41 :Definition:
43 An instance :math:`\Theta` of this class defines a perturbation parameter
45 .. math::
47 \theta \in \Theta = \operatorname{conv}(S)
49 where :math:`S \subset \mathbb{R}^{m \times n}` is a finite set of
50 *scenarios* and :math:`\operatorname{conv}` denotes the convex hull.
52 Usually, the scenarios are observed or projected realizations of the
53 uncertain data and :math:`\theta` is used to represent the data directly.
55 :Example:
57 >>> from picos.uncertain import ScenarioPerturbationSet
58 >>> scenarios = [[1, -1], [1, 1], [-1, -1], [-1, 1], [0, 0]]
59 >>> S = ScenarioPerturbationSet("s", scenarios, False); S
60 <2×1 Scenario Perturbation Set: conv({5 2×1 scenarios})>
61 >>> s = S.parameter; s
62 <2×1 Perturbation: s>
63 >>> # Compute largest sum of entries over all points in S.
64 >>> value, realization = S.worst_case(s[0] + s[1], "max")
65 >>> round(value, 4)
66 2.0
67 >>> print(realization)
68 [ 1.00e+00]
69 [ 1.00e+00]
70 <BLANKLINE>
71 """
73 def __init__(self, parameter_name, scenarios, compute_hull=None):
74 """Create a :class:`ScenarioPerturbationSet`.
76 :param str parameter_name:
77 Name of the parameter that lives in the set.
79 :param scenarios:
80 A collection of data points of same shape representing :math:`S`.
81 :type scenarios:
82 anything recognized by :class:`picos.Samples`
84 :param bool compute_hull:
85 Whether to use SciPy to compute the convex hull of the data points
86 and discard points in the interior. This can speed up the solution
87 process significantly, in particular when the scenarios come from
88 observations and when the data is low-dimensional. On the other
89 hand, when the given scenarios are known to be on the boundary of
90 their convex hull, then disabling this speeds up initialization of
91 the perturbation set. The default value of :obj:`None` means
92 :obj:`True` when SciPy is available and :obj:`False` otherwise.
93 """
94 S = Samples(scenarios)
96 if compute_hull is not False:
97 if S.dim == 1:
98 # scipy.spatial.ConvexHull does not work on scalars.
99 S = Samples([min(S._cvxopt_matrix), max(S._cvxopt_matrix)])
100 else:
101 try:
102 from scipy.spatial import qhull
103 except ModuleNotFoundError as error:
104 if compute_hull:
105 raise RuntimeError("PICOS requires SciPy to compute a "
106 "convex hull.") from error
107 else:
108 try:
109 hull = qhull.ConvexHull(cvx2np(S._cvxopt_matrix.T))
110 except qhull.QhullError as error:
111 if compute_hull:
112 raise RuntimeError("SciPy failed to compute a "
113 "convex hull.") from error
114 else:
115 S = S.select(hull.vertices.tolist())
117 # Define a convex combination of the scenarios for anticipation.
118 t = RealVariable("t", S.num, lower=0)
119 A = (S.matrix*t).reshaped(S.original_shape)
121 self._scenarios = S
122 self._convex_combination = t, A
123 self._parameter = Perturbation(self, parameter_name, S.original_shape)
125 @property
126 def scenarios(self):
127 """The registered scenarios as a :class:`~.samples.Samples` object."""
128 return self._scenarios
130 Subtype = namedtuple("Subtype", ("param_dim", "scenario_count"))
132 def _subtype(self):
133 return self.Subtype(self._parameter.dim, self._scenarios.num)
135 def __str__(self):
136 return "conv({})".format(glyphs.set("{} {} scenarios".format(
137 self._scenarios.num, glyphs.shape(self.parameter.shape))))
139 @classmethod
140 def _get_type_string_base(cls):
141 return "Scenario Perturbation Set"
143 def __repr__(self):
144 return glyphs.repr2("{} {}".format(glyphs.shape(self._parameter.shape),
145 self._get_type_string_base()), self.__str__())
147 @property
148 def distributional(self):
149 """Implement for :class:`~.perturbation.PerturbationUniverse`."""
150 return False
152 @property
153 def parameter(self):
154 """Implement for :class:`~.perturbation.PerturbationUniverse`."""
155 return self._parameter
157 def worst_case(self, scalar, direction):
158 """Implement for :class:`~.perturbation.PerturbationUniverse`."""
159 from ...modeling import Problem, SolutionFailure
161 self._check_worst_case_argument_scalar(scalar)
162 self._check_worst_case_argument_direction(direction)
164 p = self._parameter
165 x = RealVariable(p.name, p.shape)
166 f = scalar.replace_mutables({p: x})
168 self._check_worst_case_f_and_x(f, x)
170 if direction == "find" \
171 or (f.convex and direction == "min") \
172 or (f.concave and direction == "max"):
173 # Use convex optimization.
174 t, A = self._convex_combination
176 P = Problem()
177 P.set_objective(direction, f)
178 P.add_constraint(x == A)
179 P.add_constraint((t | 1) == 1)
181 try:
182 P.solve(**settings.INTERNAL_OPTIONS)
183 return f.safe_value, x.safe_value
184 except (SolutionFailure, NotValued) as error:
185 raise RuntimeError(
186 "Failed to compute {}({}) for {}.".format(direction,
187 f.string, glyphs.element(x.string, self))) from error
188 elif (f.convex and direction == "max") \
189 or (f.concave and direction == "min"):
190 # Use enumeration since at least one scenario is optimal.
191 minimizing = direction == "min"
193 value = float("inf") if minimizing else float("-inf")
194 realization = None
196 for scenario in self._scenarios._cvxopt_vectors:
197 x.value = scenario
198 f_value = f.value
200 if (minimizing and f_value < value) \
201 or (not minimizing and f_value > value):
202 value = f_value
203 realization = scenario
205 assert realization is not None
207 if realization.size == (1, 1):
208 realization = realization[0]
209 else:
210 realization = cvxopt.matrix(realization, p.shape)
212 return value, realization
213 else:
214 assert not f.convex and not f.concave
216 raise IntractableWorstCase("PICOS does not know how to compute "
217 "{}({}) for {} as the objective function seems to be neither "
218 "convex nor concave.".format(direction, scalar.string,
219 glyphs.element(p.name, self)))
222# --------------------------------------
223__all__ = api_end(_API_START, globals())