Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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# ------------------------------------------------------------------------------ 

18 

19"""Implements :class:`ScenarioPerturbationSet`.""" 

20 

21from collections import namedtuple 

22 

23import cvxopt 

24 

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 

33 

34_API_START = api_start(globals()) 

35# ------------------------------- 

36 

37 

38class ScenarioPerturbationSet(PerturbationUniverse): 

39 r"""A scenario description of a :class:`~.perturbation.Perturbation`. 

40 

41 :Definition: 

42 

43 An instance :math:`\Theta` of this class defines a perturbation parameter 

44 

45 .. math:: 

46 

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

48 

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. 

51 

52 Usually, the scenarios are observed or projected realizations of the 

53 uncertain data and :math:`\theta` is used to represent the data directly. 

54 

55 :Example: 

56 

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 """ 

72 

73 def __init__(self, parameter_name, scenarios, compute_hull=None): 

74 """Create a :class:`ScenarioPerturbationSet`. 

75 

76 :param str parameter_name: 

77 Name of the parameter that lives in the set. 

78 

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` 

83 

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) 

95 

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 ConvexHull 

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 hull = ConvexHull(cvx2np(S._cvxopt_matrix.T)) 

109 S = S.select(hull.vertices.tolist()) 

110 

111 # Define a convex combination of the scenarios for anticipation. 

112 t = RealVariable("t", S.num, lower=0) 

113 A = (S.matrix*t).reshaped(S.original_shape) 

114 

115 self._scenarios = S 

116 self._convex_combination = t, A 

117 self._parameter = Perturbation(self, parameter_name, S.original_shape) 

118 

119 @property 

120 def scenarios(self): 

121 """The registered scenarios as a :class:`~.samples.Samples` object.""" 

122 return self._scenarios 

123 

124 Subtype = namedtuple("Subtype", ("param_dim", "scenario_count")) 

125 

126 def _subtype(self): 

127 return self.Subtype(self._parameter.dim, self._scenarios.num) 

128 

129 def __str__(self): 

130 return "conv({})".format(glyphs.set("{} {} scenarios".format( 

131 self._scenarios.num, glyphs.shape(self.parameter.shape)))) 

132 

133 @classmethod 

134 def _get_type_string_base(cls): 

135 return "Scenario Perturbation Set" 

136 

137 def __repr__(self): 

138 return glyphs.repr2("{} {}".format(glyphs.shape(self._parameter.shape), 

139 self._get_type_string_base()), self.__str__()) 

140 

141 @property 

142 def distributional(self): 

143 """Implement for :class:`~.perturbation.PerturbationUniverse`.""" 

144 return False 

145 

146 @property 

147 def parameter(self): 

148 """Implement for :class:`~.perturbation.PerturbationUniverse`.""" 

149 return self._parameter 

150 

151 def worst_case(self, scalar, direction): 

152 """Implement for :class:`~.perturbation.PerturbationUniverse`.""" 

153 from ...modeling import Problem, SolutionFailure 

154 

155 self._check_worst_case_argument_scalar(scalar) 

156 self._check_worst_case_argument_direction(direction) 

157 

158 p = self._parameter 

159 x = RealVariable(p.name, p.shape) 

160 f = scalar.replace_mutables({p: x}) 

161 

162 self._check_worst_case_f_and_x(f, x) 

163 

164 if direction == "find" \ 

165 or (f.convex and direction == "min") \ 

166 or (f.concave and direction == "max"): 

167 # Use convex optimization. 

168 t, A = self._convex_combination 

169 

170 P = Problem() 

171 P.set_objective(direction, f) 

172 P.add_constraint(x == A) 

173 P.add_constraint((t | 1) == 1) 

174 

175 try: 

176 P.solve(**settings.INTERNAL_OPTIONS) 

177 return f.safe_value, x.safe_value 

178 except (SolutionFailure, NotValued) as error: 

179 raise RuntimeError( 

180 "Failed to compute {}({}) for {}.".format(direction, 

181 f.string, glyphs.element(x.string, self))) from error 

182 elif (f.convex and direction == "max") \ 

183 or (f.concave and direction == "min"): 

184 # Use enumeration since at least one scenario is optimal. 

185 minimizing = direction == "min" 

186 

187 value = float("inf") if minimizing else float("-inf") 

188 realization = None 

189 

190 for scenario in self._scenarios._cvxopt_vectors: 

191 x.value = scenario 

192 f_value = f.value 

193 

194 if (minimizing and f_value < value) \ 

195 or (not minimizing and f_value > value): 

196 value = f_value 

197 realization = scenario 

198 

199 assert realization is not None 

200 

201 if realization.size == (1, 1): 

202 realization = realization[0] 

203 else: 

204 realization = cvxopt.matrix(realization, p.shape) 

205 

206 return value, realization 

207 else: 

208 assert not f.convex and not f.concave 

209 

210 raise IntractableWorstCase("PICOS does not know how to compute " 

211 "{}({}) for {} as the objective function seems to be neither " 

212 "convex nor concave.".format(direction, scalar.string, 

213 glyphs.element(p.name, self))) 

214 

215 

216# -------------------------------------- 

217__all__ = api_end(_API_START, globals())