Coverage for picos/expressions/uncertain/pert

1# ------------------------------------------------------------------------------

4# This file is part of PICOS.

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

21from collections import namedtuple

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

32_API_START = api_start(globals())

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

36class MomentAmbiguitySet(PerturbationUniverse):

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

39 :Model of uncertainty:

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

42 class

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

50 .. math::

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

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

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

57 :Definition:

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

61 .. math::

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

73 where

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.

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.

91 .. note::

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.

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

Coverage for picos/expressions/uncertain/pert_moment.py: 75.00%

96 statements