Coverage for picos/constraints/uncertain/ucon_mom_sqnorm.py: 67.67%
133 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
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# ------------------------------------------------------------------------------
19"""Implements :class:`MomentAmbiguousSquaredNormConstraint`."""
21from collections import namedtuple
23from ... import glyphs
24from ...apidoc import api_end, api_start
25from ..constraint import Constraint, ConstraintConversion
27_API_START = api_start(globals())
28# -------------------------------
31class MomentAmbiguousSquaredNormConstraint(Constraint):
32 """A bound on a moment-ambiguous expected value of a squared norm."""
34 class DistributionallyRobustConversion(ConstraintConversion):
35 """Distributionally robust counterpart conversion."""
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
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
49 if bounded_mean or bounded_covariance:
50 yield ("var", RealVariable.make_var_type(dim=1, bnd=0), 1) # r
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
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)
62 if bounded_support:
63 # l >= 0
64 yield ("var", RealVariable.make_var_type(dim=1, bnd=1), 1)
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."
82 yield ("con", LMIConstraint.make_type(diag=(k + m + 1)), 1)
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
97 # Load the uncertain suqared norm.
98 a = con.sqnorm.x
99 B, b = a.factor_out(a.perturbation)
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
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
115 # Load the upper bound.
116 omega = con.ub
118 problem = Problem()
120 if bounded_mean or bounded_covariance:
121 r = RealVariable("__r")
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))
130 if bounded_covariance:
131 Q = SymmetricVariable("__Q", m)
132 problem.add_constraint(Q >> 0)
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
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
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
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
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
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."
187 problem.add_constraint(
188 block([["I", B, b],
189 [B.T, U, V],
190 [b.T, V.T, W]]) >> 0
191 )
193 return problem
195 def __init__(self, sqnorm, upper_bound):
196 """Construct a :class:`MomentAmbiguousSquaredNormConstraint`.
198 :param ~picos.expressions.UncertainSquaredNorm sqnorm:
199 Uncertain squared norm to upper bound the expectation of.
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
207 assert isinstance(sqnorm, UncertainSquaredNorm)
208 assert isinstance(sqnorm.universe, MomentAmbiguitySet)
209 assert isinstance(upper_bound, AffineExpression)
210 assert upper_bound.scalar
212 self.sqnorm = sqnorm
213 self.ub = upper_bound
215 super(MomentAmbiguousSquaredNormConstraint, self).__init__(
216 "Moment-ambiguous Expected Squared Norm", printSize=True)
218 Subtype = namedtuple("Subtype", ("sqnorm_argdim", "universe_subtype"))
220 def _subtype(self):
221 return self.Subtype(len(self.sqnorm.x), self.sqnorm.universe.subtype)
223 @classmethod
224 def _cost(cls, subtype):
225 return float("inf")
227 def _expression_names(self):
228 yield "sqnorm"
229 yield "ub"
231 def _str(self):
232 return glyphs.le(self.sqnorm.worst_case_string("max"), self.ub.string)
234 def _get_size(self):
235 return (1, 1)
237 def _get_slack(self):
238 return self.ub.value - self.sqnorm.worst_case_value(direction="max")
241# --------------------------------------
242__all__ = api_end(_API_START, globals())