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

180 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-03-26 07:46 +0000

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

20 

21from collections import namedtuple 

22 

23import cvxopt 

24 

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 

40 

41_API_START = api_start(globals()) 

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

43 

44 

45class ConicPerturbationSet(PerturbationUniverse): 

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

47 

48 :Definition: 

49 

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

51 

52 .. math:: 

53 

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

57 

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

63 

64 :Usage: 

65 

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

67 

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. 

79 

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

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

82 

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

86 

87 :Example: 

88 

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