Coverage for picos/expressions/uncertain/pert_scenario.py: 88.89%

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 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()) 

116 

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) 

120 

121 self._scenarios = S 

122 self._convex_combination = t, A 

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

124 

125 @property 

126 def scenarios(self): 

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

128 return self._scenarios 

129 

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

131 

132 def _subtype(self): 

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

134 

135 def __str__(self): 

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

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

138 

139 @classmethod 

140 def _get_type_string_base(cls): 

141 return "Scenario Perturbation Set" 

142 

143 def __repr__(self): 

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

145 self._get_type_string_base()), self.__str__()) 

146 

147 @property 

148 def distributional(self): 

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

150 return False 

151 

152 @property 

153 def parameter(self): 

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

155 return self._parameter 

156 

157 def worst_case(self, scalar, direction): 

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

159 from ...modeling import Problem, SolutionFailure 

160 

161 self._check_worst_case_argument_scalar(scalar) 

162 self._check_worst_case_argument_direction(direction) 

163 

164 p = self._parameter 

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

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

167 

168 self._check_worst_case_f_and_x(f, x) 

169 

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 

175 

176 P = Problem() 

177 P.set_objective(direction, f) 

178 P.add_constraint(x == A) 

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

180 

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" 

192 

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

194 realization = None 

195 

196 for scenario in self._scenarios._cvxopt_vectors: 

197 x.value = scenario 

198 f_value = f.value 

199 

200 if (minimizing and f_value < value) \ 

201 or (not minimizing and f_value > value): 

202 value = f_value 

203 realization = scenario 

204 

205 assert realization is not None 

206 

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

208 realization = realization[0] 

209 else: 

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

211 

212 return value, realization 

213 else: 

214 assert not f.convex and not f.concave 

215 

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

220 

221 

222# -------------------------------------- 

223__all__ = api_end(_API_START, globals())