Coverage for picos/constraints/uncertain/ucon_ws_sqnorm.py: 96.10%

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

20 

21from collections import namedtuple 

22 

23from ... import glyphs 

24from ...apidoc import api_end, api_start 

25from ..constraint import Constraint, ConstraintConversion 

26 

27_API_START = api_start(globals()) 

28# ------------------------------- 

29 

30 

31class WassersteinAmbiguousSquaredNormConstraint(Constraint): 

32 """A bound on a Wasserstein-ambiguous expected value of a squared norm.""" 

33 

34 class DistributionallyRobustConversion(ConstraintConversion): 

35 """Distributionally robust counterpart conversion.""" 

36 

37 @classmethod 

38 def predict(cls, subtype, options): 

39 """Implement :meth:`~.constraint.ConstraintConversion.predict`.""" 

40 from ...expressions import RealVariable, SymmetricVariable 

41 from .. import AffineConstraint, LMIConstraint 

42 

43 k = subtype.sqnorm_argdim 

44 m = subtype.universe_subtype.sample_dim 

45 N = subtype.universe_subtype.sample_num 

46 

47 yield ("var", RealVariable.make_var_type(dim=1, bnd=1), 1) # gamma 

48 yield ("var", RealVariable.make_var_type(dim=N, bnd=0), 1) # s 

49 yield ("var", SymmetricVariable.make_var_type( # U 

50 dim=(m * (m + 1)) // 2, bnd=0), 1) 

51 yield ("var", RealVariable.make_var_type(dim=m, bnd=0), 1) # u 

52 yield ("var", RealVariable.make_var_type(dim=1, bnd=0), 1) # mu 

53 

54 yield ("con", AffineConstraint.make_type(dim=1, eq=False), 1) 

55 yield ("con", LMIConstraint.make_type(diag=(k + m + 1)), 1) 

56 yield ("con", LMIConstraint.make_type(diag=(m + 1)), N) 

57 

58 @classmethod 

59 def convert(cls, con, options): 

60 """Implement :meth:`~.constraint.ConstraintConversion.convert`.""" 

61 # The recipe is found in "Robust conic optimization in 

62 # Python" (Stahlberg 2020) and extends a result in "Wasserstein 

63 # distributionally robust optimization: Theory and applications in 

64 # machine learning" (Kuhn, Esfahani, Nguyen and Shafieezadeh-Abadeh 

65 # 2019). 

66 from ...expressions import RealVariable, SymmetricVariable 

67 from ...expressions.algebra import block 

68 from ...modeling import Problem 

69 

70 problem = Problem() 

71 

72 # Load the uncertain suqared norm. 

73 a = con.sqnorm.x 

74 B, b = a.factor_out(a.perturbation) 

75 

76 # Load the ambiguity set. 

77 WAS = con.sqnorm.universe 

78 S = WAS.samples 

79 m = S.dim 

80 N = S.num 

81 w = WAS.weights 

82 eps = WAS.eps 

83 

84 # Load the upper bound. 

85 omega = con.ub 

86 

87 # Introduce auxiliary variables. 

88 gamma = RealVariable("__gamma", lower=0) 

89 s = RealVariable("__s", N) 

90 U = SymmetricVariable("__U", m) 

91 u = RealVariable("__u", m) 

92 mu = RealVariable("__mu") 

93 

94 # Compute redundant terms that appear in constraints. 

95 h1 = gamma.dupdiag(m) - U 

96 h2 = tuple(gamma*S[i] + u for i in range(N)) 

97 h3 = tuple(gamma*abs(S[i])**2 + s[i] - mu for i in range(N)) 

98 

99 # Add constraints. 

100 problem.add_constraint( 

101 gamma*eps**2 + w.T*s <= omega) 

102 problem.add_constraint( 

103 block([["I", B, b], 

104 [B.T, U, u], 

105 [b.T, u.T, mu]]) >> 0) 

106 problem.add_list_of_constraints([ 

107 block([[h1, h2[i]], 

108 [h2[i].T, h3[i]]]) >> 0 

109 for i in range(N)]) 

110 

111 return problem 

112 

113 def __init__(self, sqnorm, upper_bound): 

114 """Construct a :class:`WassersteinAmbiguousSquaredNormConstraint`. 

115 

116 :param ~picos.expressions.UncertainSquaredNorm sqnorm: 

117 Uncertain squared norm to upper bound the expectation of. 

118 

119 :param ~picos.expressions.AffineExpression upper_bound: 

120 Upper bound on the expected value. 

121 """ 

122 from ...expressions import AffineExpression, UncertainSquaredNorm 

123 from ...expressions.uncertain.pert_wasserstein import ( 

124 WassersteinAmbiguitySet) 

125 

126 assert isinstance(sqnorm, UncertainSquaredNorm) 

127 assert isinstance(sqnorm.universe, WassersteinAmbiguitySet) 

128 assert sqnorm.universe.p == 2 

129 assert isinstance(upper_bound, AffineExpression) 

130 assert upper_bound.scalar 

131 

132 self.sqnorm = sqnorm 

133 self.ub = upper_bound 

134 

135 super(WassersteinAmbiguousSquaredNormConstraint, self).__init__( 

136 "Wasserstein-ambiguous Expected Squared Norm", printSize=True) 

137 

138 Subtype = namedtuple("Subtype", ("sqnorm_argdim", "universe_subtype")) 

139 

140 def _subtype(self): 

141 return self.Subtype(len(self.sqnorm.x), self.sqnorm.universe.subtype) 

142 

143 @classmethod 

144 def _cost(cls, subtype): 

145 return float("inf") 

146 

147 def _expression_names(self): 

148 yield "sqnorm" 

149 yield "ub" 

150 

151 def _str(self): 

152 return glyphs.le(self.sqnorm.worst_case_string("max"), self.ub.string) 

153 

154 def _get_size(self): 

155 return (1, 1) 

156 

157 def _get_slack(self): 

158 return self.ub.value - self.sqnorm.worst_case_value(direction="max") 

159 

160 

161# -------------------------------------- 

162__all__ = api_end(_API_START, globals())