Coverage for picos/constraints/uncertain/ucon_scen_conic.py: 94.12%

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

20 

21import operator 

22from collections import namedtuple 

23 

24from ... import glyphs 

25from ...apidoc import api_end, api_start 

26from ..constraint import Constraint, ConstraintConversion 

27 

28_API_START = api_start(globals()) 

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

30 

31 

32class ScenarioUncertainConicConstraint(Constraint): 

33 """Conic constraint with scenario uncertainty.""" 

34 

35 class RobustConversion(ConstraintConversion): 

36 """Robust counterpart conversion.""" 

37 

38 @classmethod 

39 def predict(cls, subtype, options): 

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

41 from ...expressions import AffineExpression 

42 

43 n = subtype.dim 

44 k = subtype.scenario_count 

45 C = subtype.cone_type 

46 

47 # HACK: It's possible that a member becomes constant for some 

48 # realization of the uncertainty but we can't predict this. 

49 # This is not a problem as long as the constraint outcome of 

50 # claiming conic membership does not depend on whether the 

51 # member is constant. 

52 a = AffineExpression.make_type( 

53 shape=(n, 1), constant=False, nonneg=False) 

54 

55 yield ("con", a.predict(operator.__lshift__, C), k) 

56 

57 @classmethod 

58 def convert(cls, con, options): 

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

60 from ...modeling import Problem 

61 

62 A, b = con.element.vec.factor_out(con.element.perturbation) 

63 

64 P = Problem() 

65 

66 for s in con.element.universe.scenarios: 

67 P.add_constraint(A*s + b << con.cone) 

68 

69 return P 

70 

71 def __init__(self, element, cone): 

72 """Construct a :class:`ScenarioUncertainConicConstraint`. 

73 

74 :param ~picos.expressions.UncertainAffineExpression element: 

75 Uncertain expression constrained to be in the cone. 

76 

77 :param ~picos.expressions.Cone cone: 

78 The cone that the uncertain expression is constrained to. 

79 """ 

80 from ...expressions import Cone, UncertainAffineExpression 

81 from ...expressions.uncertain.pert_scenario import ( 

82 ScenarioPerturbationSet) 

83 

84 assert isinstance(element, UncertainAffineExpression) 

85 assert isinstance(element.universe, ScenarioPerturbationSet) 

86 assert isinstance(cone, Cone) 

87 assert cone.dim is None or len(element) == cone.dim 

88 

89 self.element = element 

90 self.cone = cone 

91 

92 super(ScenarioUncertainConicConstraint, self).__init__( 

93 "Scenario Uncertain Conic", printSize=True) 

94 

95 Subtype = namedtuple("Subtype", ("dim", "scenario_count", "cone_type")) 

96 

97 def _subtype(self): 

98 return self.Subtype( 

99 dim=len(self.element), 

100 scenario_count=self.element.universe.scenarios.num, 

101 cone_type=self.cone.type) 

102 

103 @classmethod 

104 def _cost(cls, subtype): 

105 return float("inf") 

106 

107 def _expression_names(self): 

108 yield "element" 

109 yield "cone" 

110 

111 def _str(self): 

112 return glyphs.forall( 

113 glyphs.element(self.element.string, self.cone.string), 

114 self.element.perturbation) 

115 

116 def _get_size(self): 

117 return self.element.shape 

118 

119 def _get_slack(self): 

120 from ...expressions import (NonnegativeOrthant, SecondOrderCone, 

121 PositiveSemidefiniteCone, RealVariable, 

122 RotatedSecondOrderCone) 

123 from ...expressions.data import cvxopt_hpsd 

124 

125 if isinstance(self.cone, NonnegativeOrthant): 

126 return self.element.worst_case_value(direction="min") 

127 elif isinstance(self.cone, SecondOrderCone): 

128 ub = self.element[0] 

129 norm = abs(self.element[1:]) 

130 

131 if ub.certain: 

132 return ub.value - norm.worst_case_value(direction="max") 

133 else: 

134 # TODO: Use convex optimization to compute the slack. 

135 raise NotImplementedError("Computing the slack of a scenario-" 

136 "uncertain second order conic constraint is not supported " 

137 "if the first element of the cone member is uncertain.") 

138 elif isinstance(self.cone, RotatedSecondOrderCone): 

139 ub1 = self.element[0] 

140 ub2 = self.element[1] 

141 sqnorm = abs(self.element[2:])**2 

142 

143 if ub1.certain and ub2.certain: 

144 ub1_value = ub1.value 

145 ub2_value = ub2.value 

146 

147 ub_value = ub1_value*ub2_value 

148 slack = ub_value - sqnorm.worst_case_value(direction="max") 

149 

150 if ub1_value < 0: 

151 slack = min(ub1_value, slack) 

152 

153 if ub2_value < 0: 

154 slack = min(ub2_value, slack) 

155 

156 return slack 

157 else: 

158 # TODO: Use convex optimization to compute the slack. 

159 raise NotImplementedError( 

160 "Computing the slack of a scenario-uncertain rotated second" 

161 " order conic constraint is not supported unless the first " 

162 "two elements of the cone member are certain.") 

163 elif isinstance(self.cone, PositiveSemidefiniteCone): 

164 # Devectorize the cone element to a symmetric matrix A and replace 

165 # its perturbation parameter with a real variable x. 

166 x = RealVariable("x", self.element.perturbation.shape) 

167 A = self.element.desvec.replace_mutables( 

168 {self.element.perturbation: x}) 

169 

170 # Find the least-slack matrix S by scenario enumeration. 

171 S = None 

172 for s in self.element.universe.scenarios._cvxopt_vectors: 

173 x.value = s 

174 if S is None or cvxopt_hpsd( 

175 S.safe_value_as_matrix - A.safe_value_as_matrix): 

176 S = ~A 

177 

178 # Vectorize the slack. 

179 return S.svec.safe_value 

180 else: 

181 # NOTE: This can be extended on a cone-by-cone basis if necessary. 

182 raise NotImplementedError("Computing the slack of a scenario-" 

183 "uncertain conic constraint is not supporeted for the cone {}." 

184 .format(self.cone.__class__.__name__)) 

185 

186 

187# -------------------------------------- 

188__all__ = api_end(_API_START, globals())