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

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 MomentAmbiguousSquaredNormConstraint(Constraint): 

32 """A bound on a moment-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, SOCConstraint 

42 

43 k = subtype.sqnorm_argdim 

44 m = subtype.universe_subtype.dim 

45 bounded_mean = subtype.universe_subtype.bounded_mean 

46 bounded_covariance = subtype.universe_subtype.bounded_covariance 

47 bounded_support = subtype.universe_subtype.bounded_support 

48 

49 if bounded_mean or bounded_covariance: 

50 yield ("var", RealVariable.make_var_type(dim=1, bnd=0), 1) # r 

51 

52 if bounded_mean: 

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

54 yield ("var", RealVariable.make_var_type(dim=1, bnd=0), 1) # t 

55 

56 if bounded_covariance: 

57 # Q >> 0 

58 yield ("var", SymmetricVariable.make_var_type( 

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

60 yield ("con", LMIConstraint.make_type(diag=m), 1) 

61 

62 if bounded_support: 

63 # l >= 0 

64 yield ("var", RealVariable.make_var_type(dim=1, bnd=1), 1) 

65 

66 if bounded_mean and bounded_covariance and bounded_support: 

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

68 yield ("con", SOCConstraint.make_type(argdim=m), 1) 

69 elif not bounded_mean and bounded_covariance and bounded_support: 

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

71 elif bounded_mean and not bounded_covariance and bounded_support: 

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

73 yield ("con", SOCConstraint.make_type(argdim=m), 1) 

74 elif bounded_mean and bounded_covariance and not bounded_support: 

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

76 yield ("con", SOCConstraint.make_type(argdim=m), 1) 

77 elif not (bounded_mean or bounded_covariance) and bounded_support: 

78 pass 

79 else: 

80 assert False, "Unexpected unboundedness pattern." 

81 

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

83 

84 @classmethod 

85 def convert(cls, con, options): 

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

87 # The conversion recipe is found in "Robust conic optimization in 

88 # Python" (Stahlberg 2020) and extends a result in "Models and 

89 # algorithms for distributionally robust least squares problems" 

90 # (Mehrotra and Zhang 2014). 

91 from ...expressions import (Constant, RealVariable, SecondOrderCone, 

92 SymmetricVariable) 

93 from ...expressions.algebra import block 

94 from ...expressions.data import cvxopt_principal_root 

95 from ...modeling import Problem 

96 

97 # Load the uncertain suqared norm. 

98 a = con.sqnorm.x 

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

100 

101 # Load the ambiguity set. 

102 MAS = con.sqnorm.universe 

103 m = MAS.dim 

104 mu = MAS.nominal_mean 

105 Sigma = MAS.nominal_covariance 

106 alpha = MAS.alpha 

107 beta = MAS.beta 

108 S = MAS.sample_space 

109 

110 # Determime boundedness pattern. 

111 bounded_mean = alpha is not None 

112 bounded_covariance = beta is not None 

113 bounded_support = S is not None 

114 

115 # Load the upper bound. 

116 omega = con.ub 

117 

118 problem = Problem() 

119 

120 if bounded_mean or bounded_covariance: 

121 r = RealVariable("__r") 

122 

123 if bounded_mean: 

124 q = RealVariable("__q", m) 

125 t = RealVariable("__t") 

126 sqrt_alpha = alpha**0.5 

127 sqrt_Sigma = Constant(glyphs.sqrt(Sigma.string), 

128 cvxopt_principal_root(Sigma.value_as_matrix)) 

129 

130 if bounded_covariance: 

131 Q = SymmetricVariable("__Q", m) 

132 problem.add_constraint(Q >> 0) 

133 

134 if bounded_support: 

135 l = RealVariable("__lambda", lower=0) 

136 inv_D = S.Ainv 

137 G = inv_D.T*inv_D 

138 d = S.c 

139 

140 if bounded_mean and bounded_covariance and bounded_support: 

141 # Default case. 

142 U = l*G + Q 

143 V = 0.5*q - l*G*d 

144 W = l*(d.T*G*d - 1) + r 

145 

146 problem.add_constraint( 

147 ((beta*Sigma + mu*mu.T) | Q) + mu.T*q + r + t <= omega) 

148 problem.add_constraint( 

149 (t // (sqrt_alpha*sqrt_Sigma*(2*Q*mu + q))) 

150 << SecondOrderCone()) 

151 elif not bounded_mean and bounded_covariance and bounded_support: 

152 # Unbounded mean. 

153 U = l*G + Q 

154 V = -Q*mu - l*G*d 

155 W = l*(d.T*G*d - 1) + r 

156 

157 problem.add_constraint( 

158 ((beta*Sigma - mu*mu.T) | Q) + r <= omega) 

159 elif bounded_mean and not bounded_covariance and bounded_support: 

160 # Unbounded covariance. 

161 U = l*G 

162 V = 0.5*q - l*G*d 

163 W = l*(d.T*G*d - 1) + r 

164 

165 problem.add_constraint(mu.T*q + r + t <= omega) 

166 problem.add_constraint( 

167 (t // (sqrt_alpha*sqrt_Sigma*q)) << SecondOrderCone()) 

168 elif bounded_mean and bounded_covariance and not bounded_support: 

169 # Unbounded support. 

170 U = Q 

171 V = 0.5*q 

172 W = r 

173 

174 problem.add_constraint( 

175 ((beta*Sigma + mu*mu.T) | Q) + mu.T*q + r + t <= omega) 

176 problem.add_constraint( 

177 (t // (sqrt_alpha*sqrt_Sigma*(2*Q*mu + q))) 

178 << SecondOrderCone()) 

179 elif not (bounded_mean or bounded_covariance) and bounded_support: 

180 # Unbounded mean and covariance. 

181 U = l*G 

182 V = -l*G*d 

183 W = l*(d.T*G*d - 1) + omega 

184 else: 

185 assert False, "Unexpected unboundedness pattern." 

186 

187 problem.add_constraint( 

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

189 [B.T, U, V], 

190 [b.T, V.T, W]]) >> 0 

191 ) 

192 

193 return problem 

194 

195 def __init__(self, sqnorm, upper_bound): 

196 """Construct a :class:`MomentAmbiguousSquaredNormConstraint`. 

197 

198 :param ~picos.expressions.UncertainSquaredNorm sqnorm: 

199 Uncertain squared norm to upper bound the expectation of. 

200 

201 :param ~picos.expressions.AffineExpression upper_bound: 

202 Upper bound on the expected value. 

203 """ 

204 from ...expressions import AffineExpression, UncertainSquaredNorm 

205 from ...expressions.uncertain.pert_moment import MomentAmbiguitySet 

206 

207 assert isinstance(sqnorm, UncertainSquaredNorm) 

208 assert isinstance(sqnorm.universe, MomentAmbiguitySet) 

209 assert isinstance(upper_bound, AffineExpression) 

210 assert upper_bound.scalar 

211 

212 self.sqnorm = sqnorm 

213 self.ub = upper_bound 

214 

215 super(MomentAmbiguousSquaredNormConstraint, self).__init__( 

216 "Moment-ambiguous Expected Squared Norm", printSize=True) 

217 

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

219 

220 def _subtype(self): 

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

222 

223 @classmethod 

224 def _cost(cls, subtype): 

225 return float("inf") 

226 

227 def _expression_names(self): 

228 yield "sqnorm" 

229 yield "ub" 

230 

231 def _str(self): 

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

233 

234 def _get_size(self): 

235 return (1, 1) 

236 

237 def _get_slack(self): 

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

239 

240 

241# -------------------------------------- 

242__all__ = api_end(_API_START, globals())