Coverage for picos/constraints/con_soc.py: 88.89%
63 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
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# ------------------------------------------------------------------------------
19"""Second order conce constraints."""
21from collections import namedtuple
23from .. import glyphs
24from ..apidoc import api_end, api_start
25from ..caching import cached_property
26from .constraint import ConicConstraint
28_API_START = api_start(globals())
29# -------------------------------
32class SOCConstraint(ConicConstraint):
33 """Second order (:math:`2`-norm, Lorentz) cone membership constraint."""
35 def __init__(self, normedExpression, upperBound, customString=None):
36 """Construct a :class:`SOCConstraint`.
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
47 assert isinstance(normedExpression, AffineExpression)
48 assert isinstance(upperBound, AffineExpression)
49 assert len(upperBound) == 1
51 # NOTE: len(normedExpression) == 1 is allowed even though this should
52 # rather be represented as an AbsoluteValueConstraint.
54 self.ne = normedExpression
55 self.ub = upperBound
57 super(SOCConstraint, self).__init__(
58 self._get_type_term(), customString, printSize=True)
60 def _get_type_term(self):
61 return "SOC"
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))
69 @cached_property
70 def unit_ball_form(self):
71 r"""The constraint in Euclidean norm unit ball membership form.
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
81 .. math::
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\}.
93 Therefor we can repose the constraint as two constraints:
95 .. math::
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}
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.
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.
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.
122 :raises ValueError:
123 If the upper bound is not constant, not positive, or if the matrix
124 :math:`A` is not invertible.
126 :Example:
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
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))
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))
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))
172 X = tuple(self.ne.mutables)[0]
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))
179 A = self.ne._linear_coefs[X]
180 b = self.ne.vec.cst
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))
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
195 y = RealVariable("__{}".format(X.name), X.dim)
197 return X, A_inverse*(c*y - b), y, Ball()
199 Subtype = namedtuple("Subtype", ("argdim",))
201 def _subtype(self):
202 return self.Subtype(len(self.ne))
204 @classmethod
205 def _cost(cls, subtype):
206 return subtype.argdim + 1
208 def _expression_names(self):
209 yield "ne"
210 yield "ub"
212 def _str(self):
213 return glyphs.le(glyphs.norm(self.ne.string), self.ub.string)
215 def _get_size(self):
216 return (len(self.ne) + 1, 1)
218 def _get_slack(self):
219 return self.ub.safe_value - abs(self.ne).safe_value
222# --------------------------------------
223__all__ = api_end(_API_START, globals())