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

21from collections import namedtuple

23import cvxopt

25from ... import glyphs, settings

26from ...apidoc import api_end, api_start

27from ...caching import cached_property

28from ...constraints import SOCConstraint

29from ..cone_product import ProductCone

30from ..cone_soc import SecondOrderCone

31from ..data import cvxopt_hcat, cvxopt_vcat

32from ..exp_affine import Constant

33from ..expression import NotValued

34from ..set_ball import Ball

35from ..variables import RealVariable

36from ..vectorizations import FullVectorization

37from .perturbation import Perturbation, PerturbationUniverse

38from .uexp_affine import UncertainAffineExpression

39from .uexpression import IntractableWorstCase

41_API_START = api_start(globals())

42# -------------------------------

45class ConicPerturbationSet(PerturbationUniverse):

46 r"""A conic description of a :class:`~.perturbation.Perturbation`.

48 :Definition:

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

52 .. math::

54 \theta \in \Theta = \{t \in \mathbb{R}^{m \times n}

55 \mid \exists u \in \mathbb{R}^l

56 : A\operatorname{vec}(t) + Bu + c \in K\}

58 where :math:`m, n \in \mathbb{Z}_{\geq 1}`,

59 :math:`l \in \mathbb{Z}_{\geq 0}`, :math:`A \in \mathbb{R}^{k \times mn}`,

60 :math:`B \in \mathbb{R}^{k \times l}`, :math:`c \in \mathbb{R}^k` and

61 :math:`K \subseteq \mathbb{R}^k` is a (product) cone for some

62 :math:`k \in \mathbb{Z}_{\geq 1}`.

64 :Usage:

66 Obtaining :math:`\theta` is done in a number of steps:

68 1. Create an instance of this class (see :meth:`__init__`).

69 2. Access :attr:`element` to obtain a regular, fresh

70 :class:`~.variables.RealVariable` representing :math:`t`.

71 3. Define :math:`\Theta` through any number of regular PICOS constraints

72 that depend only on :math:`t` and that have a conic representation by

73 passing the constraints to :meth:`bound`.

74 4. Call :meth:`compile` to obtain a handle to the

75 :class:`~.perturbation.Perturbation` :math:`\theta`.

76 5. You can now use :math:`\theta` to build instances of

77 :class:`~.uexp_affine.UncertainAffineExpression` and derived constraint

78 types.

80 It is best practice to assign :math:`t` to a Python variable and overwrite

81 that variable with :math:`\theta` on compilation.

83 Alternatively, you can obtain a compiled :math:`\Theta` from the factory

84 method :meth:`from_constraints` and access :math:`\theta` via

85 :attr:`parameter`.

87 :Example:

89 >>> from picos import Constant, Norm, RealVariable

90 >>> from picos.uncertain import ConicPerturbationSet

91 >>> S = ConicPerturbationSet("P", (4, 4))

92 >>> P = S.element; P # Obtain a temporary parameter to describe S with.

93 <4×4 Real Variable: P>

94 >>> S.bound(Norm(P, float("inf")) <= 1) # Confine each element to [-1,1].

95 >>> S.bound(Norm(P, 1) <= 4); S # Allow a perturbation budget of 4.

96 <4×4 Conic Perturbation Set: {P : ‖P‖_max ≤ 1 ∧ ‖P‖_sum ≤ 4}>

97 >>> P = S.compile(); P # Compile the set and obtain the actual parameter.

98 <4×4 Perturbation: P>

99 >>> A = Constant("A", range(16), (4, 4))

100 >>> U = A + P; U # Define an uncertain data matrix.

101 <4×4 Uncertain Affine Expression: A + P>

102 >>> x = RealVariable("x", 4)

103 >>> U*x # Define an uncertain affine expression.

104 <4×1 Uncertain Affine Expression: (A + P)·x>

105 """

106

107 def __init__(self, parameter_name, shape=(1, 1)):

108 """Create a :class:`ConicPerturbationSet`.

109

110 :param str parameter_name:

111 Name of the parameter that lives in the set.

112

113 :param shape:

114 The shape of a vector or matrix perturbation.

115 :type shape:

116 int or tuple or list

117 """

118 from ...modeling import Problem

119

120 self._compiled = False

121 self._element = RealVariable(parameter_name, shape)

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

123 self._bounds = Problem()

124

125 Subtype = namedtuple("Subtype", (

126 "param_dim", "cone_type", "dual_cone_type", "has_B"))

127

128 def _subtype(self):

129 return self.Subtype(

130 self._parameter.dim,

131 self.K.type,

132 self.K.dual_cone.type,

133 self.B is not None)

134

135 @classmethod

136 def from_constraints(cls, parameter_name, *constraints):

137 """Create a :class:`ConicPerturbationSet` from constraints.

138

139 The constraints must concern a single regular decision variable that

140 plays the role of the :attr:`element` :math:`t`. This variable is not

141 stored or modified and can be reused in a different context.

142

143 :param str parameter_name:

144 Name of the parameter that lives in the set.

145

146 :param constraints:

147 A parameter sequence of constraints that concern a single regular

148 decision variable whose internal vectorization is trivial (its

149 dimension must match the product over its shape) and that have a

150 conic representation.

151

152 :raises ValueError:

153 If the constraints do not all concern the same single variable.

154

155 :raises TypeError:

156 If the variable uses a nontrivial vectorization format or if the

157 constraints do not all have a conic representation.

158

159 :Example:

160

161 >>> from picos.expressions.uncertain import ConicPerturbationSet

162 >>> from picos import RealVariable

163 >>> x = RealVariable("x", 4)

164 >>> T = ConicPerturbationSet.from_constraints("t", abs(x) <= 2, x >= 0)

165 >>> print(T)

166 {t : ‖t‖ ≤ 2 ∧ t ≥ 0}

167 >>> print(repr(T.parameter))

168 <4×1 Perturbation: t>

169 """

170 T, t, seen_variable = None, None, None

171

172 for constraint in constraints:

173 if len(constraint.variables) != 1:

174 raise ValueError("The constraint {} does not concern exactly "

175 "one variable.".format(constraint))

176

177 variable = next(iter(constraint.variables))

178

179 if not isinstance(variable._vec, FullVectorization):

180 raise TypeError("The variable {} cannot be used to construct a "

181 "{} from constraints because it uses a nontrivial "

182 "vectorization format.".format(variable, cls.name))

183

184 if not seen_variable:

185 seen_variable = variable

186 T = cls(parameter_name, variable.shape)

187 t = T.element

188 elif variable is not seen_variable:

189 raise ValueError("The constraints do not concern the same "

190 "single variable (found {} and {})."

191 .format(seen_variable.name, variable.name))

192

193 T.bound(constraint.replace_mutables({variable: t}))

194

195 T.compile()

196 return T

197

198 def __str__(self):

199 return glyphs.set(glyphs.sep(self._parameter.name,

200 glyphs.and_("", "").join(str(con)

201 for con in self._bounds.constraints.values())))

202

203 @classmethod

204 def _get_type_string_base(cls):

205 return "Conic Perturbation Set"

206

207 def __repr__(self):

208 return glyphs.repr2("{} {}".format(glyphs.shape(self._element.shape),

209 self._get_type_string_base()), self.__str__())

210

211 def _forbid_compiled(self):

212 if self._compiled:

213 raise RuntimeError(

214 "The perturbation set has already been compiled.")

215

216 def _require_compiled(self):

217 if not self._compiled:

218 raise RuntimeError(

219 "The perturbation set has not yet been compiled.")

220

221 @property

222 def element(self):

223 r"""The perturbation element :math:`t` describing the set.

224

225 This is a regular :class:`~.variables.RealVariable` that you can use to

226 create constraints to pass to :meth:`bound`. You can then obtain the

227 "actual" perturbation parameter :math:`\theta` to use in expressions

228 alongside your decision variaiables using :meth:`compile`.

229

230 .. warning::

231

232 If you use this object instead of :attr:`parameter` to define a

233 decision problem then it will act as a regular decision variable,

234 which is probably not what you want.

235

236 :raises RuntimeError:

237 If the set was already compiled.

238 """

239 self._forbid_compiled()

240 return self._element

241

242 def bound(self, constraint):

243 r"""Add a constraint that bounds :math:`t`.

244

245 The constraints do not need to be conic but they need to have a *conic

246 representation*, which may involve any number of auxiliary variables.

247 For instance, given a constant *uncertainty budget* :math:`b`, you may

248 add the bound :math:`\lVert t \rVert_1 \leq b` (via

249 ``picos.Norm(t, 1) <= b``) which can be represented in conic form as

250

251 .. math::

252

253 &\exists v \in \mathbb{R}^{\operatorname{dim}(t)}

254 : -t \leq v \land t \leq v \land \mathbf{1}^Tv \leq b \\

255 \Longleftrightarrow~

256 &\exists v \in \mathbb{R}^{\operatorname{dim}(t)} :

257 \begin{pmatrix}

258 v + t \\

259 v - t \\

260 b - \mathbf{1}^Tv

261 \end{pmatrix} \in \mathbb{R}_{\geq 0}^{2\operatorname{dim}(t) + 1}.

262

263 The auxiliary variable :math:`v` then plays the role of (a slice of)

264 :math:`u` in the formal definition of :math:`\Theta`.

265

266 When you are done adding bounds, you can obtain :math:`\theta` using

267 :meth:`compile`.

268

269 :raises RuntimeError:

270 If the set was already compiled.

271 """

272 self._forbid_compiled()

273

274 vars = constraint.variables

275

276 if len(vars) != 1 or next(iter(vars)) != self._element:

277 raise ValueError(

278 "The constraint {} does not bound (only) the perturbation "

279 "element {}.".format(constraint, self._element.string))

280

281 self._bounds.add_constraint(constraint)

282

283 def compile(self, validate_feasibility=False):

284 r"""Compile the set and return :math:`\theta`.

285

286 Internally, this computes the matrices :math:`A` and :math:`B`, the

287 vector :math:`c` and the (product) cone :math:`K`.

288

289 :param bool validate_feasibility:

290 Whether to solve the feasibility problem associated with the bounds

291 on :math:`t` to verify that :math:`\Theta` is nonempty.

292

293 :returns:

294 An instance of :class:`~.perturbation.Perturbation`.

295

296 :raises RuntimeError:

297 If the set was already compiled.

298

299 :raises TypeError:

300 If the bound constraints could not be put into conic form.

301

302 :raises ValueError:

303 If :math:`\Theta` could not be verified to be nonempty (needs

304 ``validate_feasibility=True``).

305 """

306 from ...constraints import DummyConstraint

307 from ...modeling import NoStrategyFound

308

309 self._forbid_compiled()

310

311 # If no bounds are given, add a DummyConstraint to mark t free.

312 if not self._bounds.constraints:

313 self._bounds.add_constraint(DummyConstraint(self._element))

314

315 # Transform all bounds into conic form.

316 try:

317 conic_bounds = self._bounds.conic_form

318 except NoStrategyFound as error:

319 raise TypeError("Could not find a conic representation for all of "

320 "the bounds on {}.".format(self._element.string)) from error

321

322 # Validate bound feasibility if requested.

323 if validate_feasibility:

324 from ...modeling.solution import SS_OPTIMAL, SS_FEASIBLE

325

326 solution = conic_bounds.solve(

327 primals=None, apply_solution=False, **settings.INTERNAL_OPTIONS)

328

329 status = solution.primalStatus

330

331 if status not in (SS_OPTIMAL, SS_FEASIBLE):

332 raise ValueError("Could not verify that the bounds on {} are "

333 "feasible: The solver {} reports a primal solution state of"

334 " {} for the associated feasibility problem.".format(

335 self._element.string, solution.solver, status))

336

337 # Form a virtual variable u from the auxiliary variables.

338 u = [var for var in conic_bounds.variables.values()

339 if var is not self._element]

340

341 # Convert auxiliary variable bounds to additional affine constraints.

342 conic_bounds.add_list_of_constraints([

343 var.bound_constraint for var in u if var.bound_constraint])

344

345 # Reformulate bounds with respect to a single product cone K.

346 A, B, c, K = [], [], [], []

347 for constraint in conic_bounds.constraints.values():

348 member, cone = constraint.conic_membership_form

349

350 if self._element in member._linear_coefs:

351 A.append(member._linear_coefs[self._element])

352 else:

353 A.append(cvxopt.spmatrix(

354 [], [], [], size=(cone.dim, self._element.dim)))

355

356 if u:

357 B.append(cvxopt_hcat([

358 member._linear_coefs[var] if var in member._linear_coefs

359 else cvxopt.spmatrix([], [], [], size=(cone.dim, var.dim))

360 for var in u]))

361

362 c.append(member._constant_coef)

363 K.append(cone)

364

365 self._A = Constant("A", cvxopt_vcat(A))

366 self._B = Constant("B", cvxopt_vcat(B)) if u else None

367 self._c = Constant("c", cvxopt_vcat(c))

368 self._K = ProductCone(*K)

369

370 # Store the bounds in conic form to speed up worst_case.

371 self._conic_bounds = conic_bounds

372

373 self._compiled = True

374 return self._parameter

375

376 @property

377 def distributional(self):

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

379 return False

380

381 @property

382 def parameter(self):

383 r"""The perturbation parameter :math:`\theta` living in the set.

384

385 This is the object returned by :meth:`compile`.

386

387 :raises RuntimeError:

388 If the set has not been compiled.

389 """

390 self._require_compiled()

391 return self._parameter

392

393 @property

394 def A(self):

395 r"""The compiled matrix :math:`A`.

396

397 :raises RuntimeError:

398 If the set has not been compiled.

399 """

400 self._require_compiled()

401 return self._A

402

403 @property

404 def B(self):

405 r"""The compiled matrix :math:`B` or :obj:`None` if :math:`l = 0`.

406

407 :raises RuntimeError:

408 If the set has not been compiled.

409 """

410 self._require_compiled()

411 return self._B

412

413 @property

414 def c(self):

415 r"""The compiled vector :math:`c`.

416

417 :raises RuntimeError:

418 If the set has not been compiled.

419 """

420 self._require_compiled()

421 return self._c

422

423 @property

424 def K(self):

425 r"""The compiled (product) cone :math:`K`.

426

427 :raises RuntimeError:

428 If the set has not been compiled.

429 """

430 self._require_compiled()

431 return self._K

432

433 def worst_case(self, scalar, direction):

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

435 from ...modeling import SolutionFailure

436

437 self._require_compiled()

438 self._check_worst_case_argument_scalar(scalar)

439 self._check_worst_case_argument_direction(direction)

440

441 p = self._parameter

442 P = self._conic_bounds.copy()

443 x = P.variables[p.name]

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

445

446 self._check_worst_case_f_and_x(f, x)

447

448 if (direction == "min" and not f.convex) \

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

450 raise IntractableWorstCase("PICOS refuses to compute {}({}) for {} "

451 "as this is a nonconvex problem.".format(direction, f.string,

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

453

454 P.set_objective(direction, f if direction != "find" else None)

455

456 try:

457 P.solve(**settings.INTERNAL_OPTIONS)

458 return f.safe_value, x.safe_value

459 except (SolutionFailure, NotValued) as error:

460 raise RuntimeError("Failed to compute {}({}) for {}.".format(

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

462

463 @cached_property

464 def unit_ball_form(self):

465 """A recipe to repose from ellipsoidal to unit norm ball uncertainty.

466

467 If the set is :attr:`ellipsoidal`, then this is a pair ``(U, M)`` where

468 ``U`` is a :class:`~.pert_conic.UnitBallPerturbationSet` and ``M`` is a

469 dictionary mapping the old :attr:`parameter` to an affine expression of

470 the new parameter that can represent the old parameter in an expression

471 (see :meth:`~.expression.Expression.replace_mutables`).

472 The mapping ``M`` is empty if and only if the perturbation set is

473 already an instance of :class:`~.pert_conic.UnitBallPerturbationSet`.

474

475 If the uncertainty set is not ellipsoidal, then this is :obj:`None`.

476

477 See also :attr:`SOCConstraint.unit_ball_form

478 <.con_soc.SOCConstraint.unit_ball_form>`.

479

480 :Example:

481

482 >>> from picos import Problem, RealVariable, sum

483 >>> from picos.uncertain import ConicPerturbationSet

484 >>> # Create a conic perturbation set and a refinement recipe.

485 >>> T = ConicPerturbationSet("t", (2, 2))

486 >>> T.bound(abs(([[1, 2], [3, 4]] ^ T.element) + 1) <= 10)

487 >>> t = T.compile()

488 >>> U, mapping = T.unit_ball_form

489 >>> print(U)

490 {t' : ‖t'‖ ≤ 1}

491 >>> print(mapping)

492 {<2×2 Perturbation: t>: <2×2 Uncertain Affine Expression: t(t')>}

493 >>> # Define and solve a conically uncertain LP.

494 >>> X = RealVariable("X", (2, 2))

495 >>> P = Problem()

496 >>> P.set_objective("max", sum(X))

497 >>> _ = P.add_constraint(X + 2*t <= 10)

498 >>> print(repr(P.parameters["t"].universe))

499 <2×2 Conic Perturbation Set: {t : ‖[2×2]⊙t + [1]‖ ≤ 10}>

500 >>> _ = P.solve(solver="cvxopt")

501 >>> print(X)

502 [-8.00e+00 1.00e+00]

503 [ 4.00e+00 5.50e+00]

504 >>> # Refine the problem to a unit ball uncertain LP.

505 >>> Q = Problem()

506 >>> Q.set_objective("max", sum(X))

507 >>> _ = Q.add_constraint(X + 2*mapping[t] <= 10)

508 >>> print(repr(Q.parameters["t'"].universe))

509 <2×2 Unit Ball Perturbation Set: {t' : ‖t'‖ ≤ 1}>

510 >>> _ = Q.solve(solver="cvxopt")

511 >>> print(X)

512 [-8.00e+00 1.00e+00]

513 [ 4.00e+00 5.50e+00]

514 """

515 self._require_compiled()

516

517 if self._B is not None:

518 return None

519

520 K = self._K.refined

521

522 if not isinstance(K, SecondOrderCone):

523 return None

524

525 C = self._A*self._element.vec + self._c << K

526

527 assert isinstance(C, SOCConstraint)

528

529 try:

530 X, aff_y, y, _ = C.unit_ball_form

531 except ValueError:

532 return None

533

534 assert X is self._element

535

536 U = UnitBallPerturbationSet("{}'".format(X.name), X.shape)

537 u = U.parameter

538 replacement = UncertainAffineExpression("{}({})".format(X.name, u.name),

539 X.shape, {(u,): aff_y._linear_coefs[y], (): aff_y._constant_coef})

540

541 return U, {self._parameter: replacement}

542

543 @property

544 def ellipsoidal(self):

545 """Whether the perturbation set is an ellipsoid.

546

547 If this is true, then a :attr:`unit_ball_form` is available.

548 """

549 return bool(self.unit_ball_form)

550

551

552class UnitBallPerturbationSet(ConicPerturbationSet):

553 r"""Represents perturbation in an Euclidean or Frobenius unit norm ball.

554

555 This is a :class:`~.pert_conic.ConicPerturbationSet` with fixed form

556

557 .. math::

558

559 \{t \in \mathbb{R}^{m \times n} \mid \lVert t \rVert_F \leq 1\}.

560

561 After initialization, you can obtain the parameter using

562 :attr:`~.pert_conic.ConicPerturbationSet.parameter`.

563 """

564

565 def __init__(self, parameter_name, shape=(1, 1)):

566 """See :meth:`ConicPerturbationSet.__init__`."""

567 ConicPerturbationSet.__init__(self, parameter_name, shape)

568 self.bound(self._element << Ball())

569 self.compile()

570

571 @classmethod

572 def _get_type_string_base(cls):

573 return "Unit Ball Perturbation Set"

574

575 @property

576 def unit_ball_form(self):

577 """Overwrite :attr:`ConicPerturbationSet.unit_ball_form`."""

578 return self, {}

579

580

581# --------------------------------------

582__all__ = api_end(_API_START, globals())

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

180 statements