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:`MomentAmbiguitySet`.""" 

20 

21from collections import namedtuple 

22 

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 

31 

32_API_START = api_start(globals()) 

33# ------------------------------- 

34 

35 

36class MomentAmbiguitySet(PerturbationUniverse): 

37 r"""A moment-uncertain description of a random perturbation parameter. 

38 

39 :Model of uncertainty: 

40 

41 As a distributional ambiguity set, an instance :math:`\mathcal{P}` of this 

42 class 

43 

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 

49 

50 .. math:: 

51 

52 \mathop{(\max\;\textit{or}\;\min)}_{\Xi \in \mathcal{P}} 

53 \mathbb{E}_\Xi[f(x, \xi)] 

54 

55 for a selection of functions :math:`f` and a decision variable :math:`x`. 

56 

57 :Definition: 

58 

59 Formally, this class can describe ambiguity sets of the form 

60 

61 .. math:: 

62 

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\} 

72 

73 where 

74 

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. 

86 

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. 

90 

91 .. note:: 

92 

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. 

99 

100 :Supported functions: 

101 

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. 

113 

114 :Example: 

115 

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. 

120 

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

150 

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

154 

155 :param str parameter_name: 

156 Name of the random parameter :math:`\xi`. 

157 

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` 

163 

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` 

169 

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` 

175 

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. 

180 

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. 

185 

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] 

202 

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 

209 

210 if not isinstance(nominal_mean, AffineExpression) \ 

211 or not nominal_mean.constant: 

212 raise TypeError("The nominal mean must be a real constant.") 

213 

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

218 

219 self._mean = nominal_mean 

220 

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 

227 

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

232 

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

237 

238 cov_value = nominal_covariance.safe_value_as_matrix 

239 

240 if not cvxopt_hpd(cov_value): 

241 raise ValueError("The nominal covariance matrix is not symmetric " 

242 "positive definite.") 

243 

244 self._cov = nominal_covariance 

245 

246 # Compute the inverse of the nominal covariance matrix. 

247 self._cov_inv = AffineExpression.from_constant( 

248 cvxopt_inverse(cov_value), shape=(n, n)) 

249 

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

257 

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

265 

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

272 

273 self._ss = None 

274 elif isinstance(sample_space, Ellipsoid): 

275 if sample_space.dim != self._dim: 

276 raise TypeError("Invalid sample space dimension.") 

277 

278 self._ss = sample_space 

279 else: 

280 raise TypeError("Invalid sample space type.") 

281 

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

286 

287 # Create the perturbation parameter. 

288 self._parameter = Perturbation(self, parameter_name, shape) 

289 

290 @property 

291 def dim(self): 

292 """The dimension :math:`n` of the sample space.""" 

293 return self._dim 

294 

295 @property 

296 def nominal_mean(self): 

297 r"""The nominal mean :math:`\mu` of the ambiguous distribution.""" 

298 return self._mean 

299 

300 @property 

301 def nominal_covariance(self): 

302 r"""The nominal covariance matrix :math:`\Sigma`.""" 

303 return self._cov 

304 

305 @property 

306 def alpha(self): 

307 r"""The parameter :math:`\alpha`. 

308 

309 A value of :obj:`None` denotes :math:`\alpha \to \infty`. 

310 """ 

311 return self._alpha 

312 

313 @property 

314 def beta(self): 

315 r"""The parameter :math:`\beta`. 

316 

317 A value of :obj:`None` denotes :math:`\beta \to \infty`. 

318 """ 

319 return self._beta 

320 

321 @property 

322 def sample_space(self): 

323 r"""The sample space (bound on support) :math:`\mathcal{S}`. 

324 

325 A value of :obj:`None` means :math:`\mathcal{S} = \mathbb{R}^n`. 

326 """ 

327 return self._ss 

328 

329 Subtype = namedtuple("Subtype", ( 

330 "dim", 

331 "bounded_mean", 

332 "bounded_covariance", 

333 "bounded_support")) 

334 

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) 

341 

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) 

349 

350 @classmethod 

351 def _get_type_string_base(cls): 

352 return "Moment Ambiguity Set" 

353 

354 def __repr__(self): 

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

356 self._get_type_string_base()), self.__str__()) 

357 

358 @property 

359 def distributional(self): 

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

361 return True 

362 

363 @property 

364 def parameter(self): 

365 r"""The random perturbation parameter :math:`\xi`.""" 

366 return self._parameter 

367 

368 

369# -------------------------------------- 

370__all__ = api_end(_API_START, globals())