Coverage for picos/constraints/uncertain/ucon_mom_pwl.py: 72.90%
155 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:`MomentAmbiguousExtremumAffineConstraint`."""
21from collections import namedtuple
23from ... import glyphs
24from ...apidoc import api_end, api_start
25from ...caching import cached_property
26from ..constraint import Constraint, ConstraintConversion
28_API_START = api_start(globals())
29# -------------------------------
32class MomentAmbiguousExtremumAffineConstraint(Constraint):
33 """A bound on a moment-ambiguous expected value of a piecewise function."""
35 class DistributionallyRobustConversion(ConstraintConversion):
36 """Distributionally robust counterpart conversion."""
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
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
50 if bounded_mean or bounded_covariance:
51 yield ("var", RealVariable.make_var_type(dim=1, bnd=0), 1) # r
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
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)
63 if bounded_support:
64 # l >= 0
65 yield ("var", RealVariable.make_var_type(dim=k, bnd=k), 1)
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."
83 yield ("con", LMIConstraint.make_type(diag=(m + 1)), k)
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
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)
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
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)
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
127 # Load the upper bound.
128 omega = con.rhs
130 problem = Problem()
132 if bounded_mean or bounded_covariance:
133 r = RealVariable("__r")
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))
142 if bounded_covariance:
143 Q = SymmetricVariable("__Q", m)
144 problem.add_constraint(Q >> 0)
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
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]
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]
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]
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]
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]
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."
215 problem.add_list_of_constraints([
216 block([[U[i], V[i]],
217 [V[i].T, W[i]]]) >> 0
218 for i in K])
220 return problem
222 def __init__(self, extremum, relation, rhs):
223 """Construct a :class:`MomentAmbiguousExtremumAffineConstraint`.
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)
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
243 self.extremum = extremum
244 self.relation = relation
245 self.rhs = rhs
247 super(MomentAmbiguousExtremumAffineConstraint, self).__init__(
248 "Moment-ambiguous Piecewise Linear Expectation", printSize=True)
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)
258 Subtype = namedtuple("Subtype", (
259 "extremum_argnum",
260 "universe_subtype"))
262 def _subtype(self):
263 return self.Subtype(
264 self.extremum.argnum,
265 self.extremum.universe.subtype)
267 @classmethod
268 def _cost(cls, subtype):
269 return float("inf")
271 def _expression_names(self):
272 yield "extremum"
273 yield "rhs"
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)
283 def _get_size(self):
284 return (1, 1)
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
293# --------------------------------------
294__all__ = api_end(_API_START, globals())