Coverage for picos/constraints/con_soc.py: 88.89%

1# ------------------------------------------------------------------------------ 

2# Copyright (C) 2018-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"""Second order conce constraints.""" 

20 

21from collections import namedtuple 

22 

23from .. import glyphs 

24from ..apidoc import api_end, api_start 

25from ..caching import cached_property 

26from .constraint import ConicConstraint 

27 

28_API_START = api_start(globals()) 

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

30 

31 

32class SOCConstraint(ConicConstraint): 

33 """Second order (:math:`2`-norm, Lorentz) cone membership constraint.""" 

34 

35 def __init__(self, normedExpression, upperBound, customString=None): 

36 """Construct a :class:`SOCConstraint`. 

37 

38 :param ~picos.expressions.AffineExpression normedExpression: 

39 Expression under the norm. 

40 :param ~picos.expressions.AffineExpression upperBound: 

41 Upper bound on the normed expression. 

42 :param str customString: 

43 Optional string description. 

44 """ 

45 from ..expressions import AffineExpression 

46 

47 assert isinstance(normedExpression, AffineExpression) 

48 assert isinstance(upperBound, AffineExpression) 

49 assert len(upperBound) == 1 

50 

51 # NOTE: len(normedExpression) == 1 is allowed even though this should 

52 # rather be represented as an AbsoluteValueConstraint. 

53 

54 self.ne = normedExpression 

55 self.ub = upperBound 

56 

57 super(SOCConstraint, self).__init__( 

58 self._get_type_term(), customString, printSize=True) 

59 

60 def _get_type_term(self): 

61 return "SOC" 

62 

63 @cached_property 

64 def conic_membership_form(self): 

65 """Implement for :class:`~.constraint.ConicConstraint`.""" 

66 from ..expressions import SecondOrderCone 

67 return (self.ub // self.ne.vec), SecondOrderCone(dim=(len(self.ne) + 1)) 

68 

69 @cached_property 

70 def unit_ball_form(self): 

71 r"""The constraint in Euclidean norm unit ball membership form. 

72 

73 If this constraint has the form :math:`\lVert E(X) \rVert_F \leq c` with 

74 :math:`c > 0` constant and :math:`E(X)` an affine expression of a single 

75 (matrix) variable :math:`X` with :math:`y := \operatorname{vec}(E(X)) = 

76 A\operatorname{vec}(X) + b` for some invertible matrix :math:`A` and a 

77 vector :math:`b`, then we have :math:`\operatorname{vec}(X) = A^{-1}(y - 

78 b)` and we can write the elementwise vectorization of the constraint's 

79 feasible region as 

80 

81 .. math:: 

82 

83 &\left\{\operatorname{vec}(X) \mid 

84 \lVert E(X) \rVert_F \leq c \right\} \\ 

85 =~&\left\{\operatorname{vec}(X) \mid 

86 \lVert A\operatorname{vec}(X) + b \rVert_2 \leq c \right\} \\ 

87 =~&\left\{A^{-1}(y - b) \mid \lVert y \rVert_2 \leq c \right\} \\ 

88 =~&\left\{ 

89 A^{-1}(y - b) \mid \lVert c^{-1}y \rVert_2 \leq 1 

90 \right\} \\ 

91 =~&\left\{A^{-1}(cy - b) \mid \lVert y \rVert_2 \leq 1 \right\}. 

92 

93 Therefor we can repose the constraint as two constraints: 

94 

95 .. math:: 

96 

97 \lVert E(X) \rVert_F \leq c 

98 \Longleftrightarrow 

99 \exists y : 

100 \begin{cases} 

101 \operatorname{vec}(X) = A^{-1}(cy - b) \\ 

102 \lVert y \rVert_2 \leq 1. 

103 \end{cases} 

104 

105 This method returns the quadruple :math:`(X, A^{-1}(cy - b), y, B)` 

106 where :math:`y` is a fresh real variable vector (the same for subsequent 

107 calls) and :math:`B` is the Euclidean norm unit ball. 

108 

109 :returns: 

110 A quadruple ``(X, aff_y, y, B)`` of type 

111 (:class:`~.variables.BaseVariable`, 

112 :class:`~.exp_affine.AffineExpression`, 

113 :class:`~.variables.RealVariable`, :class:`~picos.Ball`) such that 

114 the two constraints ``X.vec == aff_y`` and ``y << B`` combined are 

115 equivalent to this one. 

116 

117 :raises NotImplementedError: 

118 If the expression under the norm does not reference exactly one 

119 variable or if that variable does not use a trivial vectorization 

120 format internally. 

121 

122 :raises ValueError: 

123 If the upper bound is not constant, not positive, or if the matrix 

124 :math:`A` is not invertible. 

125 

126 :Example: 

127 

128 >>> import picos 

129 >>> A = picos.Constant("A", [[2, 0], 

130 ... [0, 1]]) 

131 >>> x = picos.RealVariable("x", 2) 

132 >>> P = picos.Problem() 

133 >>> P.set_objective("max", picos.sum(x)) 

134 >>> C = P.add_constraint(abs(A*x + 1) <= 10) 

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

136 >>> print(x) 

137 [ 1.74e+00] 

138 [ 7.94e+00] 

139 >>> Q = picos.Problem() 

140 >>> Q.set_objective("max", picos.sum(x)) 

141 >>> x, aff_y, y, B, = C.unit_ball_form 

142 >>> _ = Q.add_constraint(x == aff_y) 

143 >>> _ = Q.add_constraint(y << B) 

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

145 >>> print(x) 

146 [ 1.74e+00] 

147 [ 7.94e+00] 

148 >>> round(abs(P.value - Q.value), 4) 

149 0.0 

150 >>> round(y[0]**2 + y[1]**2, 4) 

151 1.0 

152 """ 

153 from ..expressions import Ball, RealVariable 

154 from ..expressions.data import cvxopt_inverse 

155 from ..expressions.vectorizations import FullVectorization 

156 

157 if len(self.ne.mutables) != 1: 

158 raise NotImplementedError("Unit ball membership form is only " 

159 "supported for second order conic constraints whose normed " 

160 "expression depends on exactly one mutable; found {} for {}." 

161 .format(len(self.ne.mutables), self)) 

162 

163 if not self.ub.constant: 

164 raise ValueError("The upper bound is not constant, so no unit ball " 

165 "form exists for {}.".format(self)) 

166 

167 c = self.ub 

168 if c.value <= 0: 

169 raise ValueError("The upper bound is not positive, so no unit ball " 

170 "form exists for {}.".format(self)) 

171 

172 X = tuple(self.ne.mutables)[0] 

173 

174 if not isinstance(X._vec, FullVectorization): 

175 raise NotImplementedError( 

176 "The variable {} does not use a trivial vectorization format, " 

177 "so no unit ball form exists for {}.".format(X.name, self)) 

178 

179 A = self.ne._linear_coefs[X] 

180 b = self.ne.vec.cst 

181 

182 if not A.size[0] == A.size[1]: 

183 raise ValueError("The dimensions dim({}) = {} and dim({}) = {} do " 

184 "not match, so no unit ball form exists for {}.".format( 

185 X.name, X.dim, self.ne.string, len(self.ne), self)) 

186 

187 try: 

188 A_inverse = cvxopt_inverse(A) 

189 except ValueError as error: 

190 raise ValueError( 

191 "The linear operator applied to {} to form the linear part of " 

192 "{} is not bijective, so no unit ball form exists for {}." 

193 .format(X.name, self.ne.string, self)) from error 

194 

195 y = RealVariable("__{}".format(X.name), X.dim) 

196 

197 return X, A_inverse*(c*y - b), y, Ball() 

198 

199 Subtype = namedtuple("Subtype", ("argdim",)) 

200 

201 def _subtype(self): 

202 return self.Subtype(len(self.ne)) 

203 

204 @classmethod 

205 def _cost(cls, subtype): 

206 return subtype.argdim + 1 

207 

208 def _expression_names(self): 

209 yield "ne" 

210 yield "ub" 

211 

212 def _str(self): 

213 return glyphs.le(glyphs.norm(self.ne.string), self.ub.string) 

214 

215 def _get_size(self): 

216 return (len(self.ne) + 1, 1) 

217 

218 def _get_slack(self): 

219 return self.ub.safe_value - abs(self.ne).safe_value 

220 

221 

222# -------------------------------------- 

223__all__ = api_end(_API_START, globals())