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

20 

21from collections import namedtuple 

22 

23from ... import glyphs 

24from ...apidoc import api_end, api_start 

25from ...caching import cached_property 

26from ..constraint import Constraint, ConstraintConversion 

27 

28_API_START = api_start(globals()) 

29# ------------------------------- 

30 

31 

32class MomentAmbiguousExtremumAffineConstraint(Constraint): 

33 """A bound on a moment-ambiguous expected value of a piecewise function.""" 

34 

35 class DistributionallyRobustConversion(ConstraintConversion): 

36 """Distributionally robust counterpart conversion.""" 

37 

38 @classmethod 

39 def predict(cls, subtype, options): 

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

41 from ...expressions import RealVariable, SymmetricVariable 

42 from .. import AffineConstraint, LMIConstraint, SOCConstraint 

43 

44 k = subtype.extremum_argnum 

45 m = subtype.universe_subtype.dim 

46 bounded_mean = subtype.universe_subtype.bounded_mean 

47 bounded_covariance = subtype.universe_subtype.bounded_covariance 

48 bounded_support = subtype.universe_subtype.bounded_support 

49 

50 if bounded_mean or bounded_covariance: 

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

52 

53 if bounded_mean: 

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

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

56 

57 if bounded_covariance: 

58 # Q >> 0 

59 yield ("var", SymmetricVariable.make_var_type( 

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

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

62 

63 if bounded_support: 

64 # l >= 0 

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

66 

67 if bounded_mean and bounded_covariance and bounded_support: 

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

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

70 elif not bounded_mean and bounded_covariance and bounded_support: 

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

72 elif bounded_mean and not bounded_covariance and bounded_support: 

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

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

75 elif bounded_mean and bounded_covariance and not bounded_support: 

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

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

78 elif not (bounded_mean or bounded_covariance) and bounded_support: 

79 pass 

80 else: 

81 assert False, "Unexpected unboundedness pattern." 

82 

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

84 

85 @classmethod 

86 def convert(cls, con, options): 

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

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

89 # Python" (Stahlberg 2020) and extends a result in "Distributionally 

90 # robust optimization under moment uncertainty with application to 

91 # data-driven problems" (Delage and Ye 2010). 

92 from ...expressions import ( 

93 Constant, RealVariable, SecondOrderCone, SymmetricVariable) 

94 from ...expressions.algebra import block 

95 from ...expressions.data import cvxopt_principal_root 

96 from ...modeling import Problem 

97 

98 # Load the constraint in maximum form. 

99 con = con.maximum_form 

100 assert isinstance(con, MomentAmbiguousExtremumAffineConstraint) 

101 assert con.relation == con.LE 

102 k = con.extremum.argnum 

103 K = range(k) 

104 

105 # Load the ambiguity set. 

106 MAS = con.extremum.universe 

107 xi = MAS.parameter 

108 m = MAS.dim 

109 mu = MAS.nominal_mean 

110 Sigma = MAS.nominal_covariance 

111 alpha = MAS.alpha 

112 beta = MAS.beta 

113 S = MAS.sample_space 

114 

115 # Load the uncertain extremum as max(h[i].T*xi + eta[i] for i in K). 

116 zero = Constant(0, shape=(1, m)) 

117 hT, eta = zip(*( 

118 (zero, a) if a.certain else a.factor_out(xi) 

119 for a in con.extremum.expressions)) 

120 h = tuple(hiT.T for hiT in hT) 

121 

122 # Determime boundedness pattern. 

123 bounded_mean = alpha is not None 

124 bounded_covariance = beta is not None 

125 bounded_support = S is not None 

126 

127 # Load the upper bound. 

128 omega = con.rhs 

129 

130 problem = Problem() 

131 

132 if bounded_mean or bounded_covariance: 

133 r = RealVariable("__r") 

134 

135 if bounded_mean: 

136 q = RealVariable("__q", m) 

137 t = RealVariable("__t") 

138 sqrt_alpha = alpha**0.5 

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

140 cvxopt_principal_root(Sigma.value_as_matrix)) 

141 

142 if bounded_covariance: 

143 Q = SymmetricVariable("__Q", m) 

144 problem.add_constraint(Q >> 0) 

145 

146 if bounded_support: 

147 l = RealVariable("__lambda", k, lower=0) 

148 l = [l[i] for i in K] # Cache element access. 

149 inv_D = S.Ainv 

150 G = inv_D.T*inv_D 

151 d = S.c 

152 Gd = G*d 

153 dTGd_minus_one = d.T*Gd - 1 

154 

155 if bounded_mean and bounded_covariance and bounded_support: 

156 # Default case. 

157 U = [l[i]*G + Q for i in K] 

158 V = [0.5*(q - h[i]) - l[i]*Gd for i in K] 

159 W = [l[i]*dTGd_minus_one - eta[i] + r for i in K] 

160 

161 problem.add_constraint( 

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

163 problem.add_constraint( 

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

165 << SecondOrderCone()) 

166 elif not bounded_mean and bounded_covariance and bounded_support: 

167 # Unbounded mean. 

168 U = [l[i]*G + Q for i in K] 

169 V = [-(Q*mu + 0.5*h[i] + l[i]*Gd) for i in K] 

170 W = [l[i]*dTGd_minus_one - eta[i] + r for i in K] 

171 

172 problem.add_constraint( 

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

174 elif bounded_mean and not bounded_covariance and bounded_support: 

175 # Unbounded covariance. 

176 U = [l[i]*G for i in K] 

177 V = [0.5*(q - h[i]) - l[i]*Gd for i in K] 

178 W = [l[i]*dTGd_minus_one - eta[i] + r for i in K] 

179 

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

181 problem.add_constraint( 

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

183 elif bounded_mean and bounded_covariance and not bounded_support: 

184 # Unbounded support. 

185 U = [Q for i in K] 

186 V = [0.5*(q - h[i]) for i in K] 

187 W = [r - eta[i] for i in K] 

188 

189 problem.add_constraint( 

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

191 problem.add_constraint( 

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

193 << SecondOrderCone()) 

194 elif not (bounded_mean or bounded_covariance) and bounded_support: 

195 # Unbounded mean and covariance. 

196 U = [l[i]*G for i in K] 

197 V = [-(0.5*h[i] + l[i]*Gd) for i in K] 

198 W = [l[i]*dTGd_minus_one - eta[i] + omega for i in K] 

199 

200 # TODO: There is also an SOCP representation (sketched below); 

201 # use it once there is sufficient test coverage to support 

202 # the special treatment. 

203 # z = [RealVariable("__z_{}".format(i), m) for i in K] 

204 # zeta = [RealVariable("__z_{}".format(i), 1) for i in K] 

205 # problem.add_list_of_constraints( 

206 # [zeta[i] - d.T*inv_D.T*z[i] + eta[i] <= omega for i in K]) 

207 # problem.add_list_of_constraints( 

208 # [inv_D.T*z[i] + h[i] == 0 for i in K]) 

209 # problem.add_list_of_constraints( 

210 # [abs(z[i]) <= zeta[i] for i in K]) 

211 # return problem 

212 else: 

213 assert False, "Unexpected unboundedness pattern." 

214 

215 problem.add_list_of_constraints([ 

216 block([[U[i], V[i]], 

217 [V[i].T, W[i]]]) >> 0 

218 for i in K]) 

219 

220 return problem 

221 

222 def __init__(self, extremum, relation, rhs): 

223 """Construct a :class:`MomentAmbiguousExtremumAffineConstraint`. 

224 

225 :param ~picos.expressions.RandomExtremumAffine extremum: 

226 Left hand side expression. 

227 :param str relation: 

228 Constraint relation symbol. 

229 :param ~picos.expressions.AffineExpression rhs: 

230 Right hand side expression. 

231 """ 

232 from ...expressions import (AffineExpression, MomentAmbiguitySet, 

233 RandomMaximumAffine, RandomMinimumAffine) 

234 

235 if relation == self.LE: 

236 assert isinstance(extremum, RandomMaximumAffine) 

237 else: 

238 assert isinstance(extremum, RandomMinimumAffine) 

239 assert isinstance(extremum.universe, MomentAmbiguitySet) 

240 assert isinstance(rhs, AffineExpression) 

241 assert rhs.scalar 

242 

243 self.extremum = extremum 

244 self.relation = relation 

245 self.rhs = rhs 

246 

247 super(MomentAmbiguousExtremumAffineConstraint, self).__init__( 

248 "Moment-ambiguous Piecewise Linear Expectation", printSize=True) 

249 

250 @cached_property 

251 def maximum_form(self): 

252 """The constraint posed as an upper bound on an expected maximum.""" 

253 if self.relation == self.LE: 

254 return self 

255 else: 

256 return self.__class__(-self.extremum, self.LE, -self.rhs) 

257 

258 Subtype = namedtuple("Subtype", ( 

259 "extremum_argnum", 

260 "universe_subtype")) 

261 

262 def _subtype(self): 

263 return self.Subtype( 

264 self.extremum.argnum, 

265 self.extremum.universe.subtype) 

266 

267 @classmethod 

268 def _cost(cls, subtype): 

269 return float("inf") 

270 

271 def _expression_names(self): 

272 yield "extremum" 

273 yield "rhs" 

274 

275 def _str(self): 

276 if self.relation == self.LE: 

277 return glyphs.le( 

278 self.extremum.worst_case_string("max"), self.rhs.string) 

279 else: 

280 return glyphs.ge( 

281 self.extremum.worst_case_string("min"), self.rhs.string) 

282 

283 def _get_size(self): 

284 return (1, 1) 

285 

286 def _get_slack(self): 

287 if self.relation == self.LE: 

288 return self.rhs.safe_value - self.extremum.worst_case_value("max") 

289 else: 

290 return self.extremum.worst_case_value("min") - self.rhs.safe_value 

291 

292 

293# -------------------------------------- 

294__all__ = api_end(_API_START, globals())