"""
Copyright 2021 the CVXPY developers
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""
from typing import List, Tuple
import numpy as np
from cvxpy.constraints.cones import Cone
from cvxpy.expressions import cvxtypes
from cvxpy.utilities import scopes
[docs]class PowCone3D(Cone):
"""
An object representing a collection of 3D power cone constraints
x[i]**alpha[i] * y[i]**(1-alpha[i]) >= |z[i]| for all i
x >= 0, y >= 0
If the parameter alpha is a scalar, it will be promoted to
a vector matching the (common) sizes of x, y, z. The numeric
value of alpha (or its components, in the vector case) must
be a number in the open interval (0, 1).
We store flattened representations of the arguments (x, y, z,
and alpha) as Expression objects. We construct dual variables
with respect to these flattened representations.
"""
def __init__(self, x, y, z, alpha, constr_id=None) -> None:
Expression = cvxtypes.expression()
self.x = Expression.cast_to_const(x)
self.y = Expression.cast_to_const(y)
self.z = Expression.cast_to_const(z)
for val in [self.x, self.y, self.z]:
if not (val.is_affine() and val.is_real()):
raise ValueError('All arguments must be affine and real.')
alpha = Expression.cast_to_const(alpha)
if alpha.is_scalar():
if self.x.shape:
alpha = cvxtypes.promote()(alpha, self.x.shape)
else:
# when `alpha` is a naked float, it has to be cast into a
# 1-D array to be compatible with downstream (vectorized)
# processing
alpha = cvxtypes.promote()(alpha, (1,))
self.alpha = alpha
if np.any(self.alpha.value <= 0) or np.any(self.alpha.value >= 1):
msg = "Argument alpha must have entries in the open interval (0, 1)."
raise ValueError(msg)
if alpha.shape == (1,):
arg_shapes = [self.x.shape, self.y.shape, self.z.shape, ()]
else:
arg_shapes = [self.x.shape, self.y.shape, self.z.shape, self.alpha.shape]
if any(arg_shapes[0] != s for s in arg_shapes[1:]):
msg = ("All arguments must have the same shapes. Provided arguments have"
"shapes %s" % str(arg_shapes))
raise ValueError(msg)
super(PowCone3D, self).__init__([self.x, self.y, self.z],
constr_id)
def __str__(self) -> str:
return "Pow3D(%s, %s, %s; %s)" % (self.x, self.y, self.z, self.alpha)
@property
def residual(self):
# TODO: The projection should be implemented directly.
from cvxpy import Minimize, Problem, Variable, hstack, norm2
if self.x.value is None or self.y.value is None or self.z.value is None:
return None
x = Variable(self.x.shape)
y = Variable(self.y.shape)
z = Variable(self.z.shape)
constr = [PowCone3D(x, y, z, self.alpha)]
obj = Minimize(norm2(hstack([x, y, z]) -
hstack([self.x.value, self.y.value, self.z.value])))
problem = Problem(obj, constr)
return problem.solve(solver='SCS', eps=1e-8)
def get_data(self):
return [self.alpha, self.id]
def is_imag(self) -> bool:
return False
def is_complex(self) -> bool:
return False
@property
def size(self) -> int:
return 3 * self.num_cones()
def num_cones(self):
return self.x.size
def cone_sizes(self) -> List[int]:
return [3]*self.num_cones()
[docs] def is_dcp(self, dpp: bool = False) -> bool:
if dpp:
with scopes.dpp_scope():
args_ok = all(arg.is_affine() for arg in self.args)
exps_ok = not isinstance(self.alpha, cvxtypes.parameter())
return args_ok and exps_ok
return all(arg.is_affine() for arg in self.args)
def is_dgp(self, dpp: bool = False) -> bool:
return False
def is_dqcp(self) -> bool:
return self.is_dcp()
@property
def shape(self) -> Tuple[int, ...]:
s = (3,) + self.x.shape
# Note: this can be a 3-tuple of x.ndim == 2.
return s
def save_dual_value(self, value) -> None:
value = np.reshape(value, newshape=(3, -1))
dv0 = np.reshape(value[0, :], newshape=self.x.shape)
dv1 = np.reshape(value[1, :], newshape=self.y.shape)
dv2 = np.reshape(value[2, :], newshape=self.z.shape)
self.dual_variables[0].save_value(dv0)
self.dual_variables[1].save_value(dv1)
self.dual_variables[2].save_value(dv2)
# TODO: figure out why the reshaping had to be done differently,
# relative to ExpCone constraints.
def _dual_cone(self, *args):
"""Implements the dual cone of PowCone3D See Pg 85
of the MOSEK modelling cookbook for more information"""
if args is None:
PowCone3D(self.dual_variables[0]/self.alpha, self.dual_variables[1]/(1-self.alpha),
self.dual_variables[2], self.alpha)
else:
# some assertions for verifying `args`
args_shapes = [arg.shape for arg in args]
instance_args_shapes = [arg.shape for arg in self.args]
assert len(args) == len(self.args)
assert args_shapes == instance_args_shapes
return PowCone3D(args[0]/self.alpha, args[1]/(1-self.alpha),
args[2], self.alpha)
[docs]class PowConeND(Cone):
"""
Represents a collection of N-dimensional power cone constraints
that is *mathematically* equivalent to the following code
snippet (which makes incorrect use of numpy functions on cvxpy
objects):
np.prod(np.power(W, alpha), axis=axis) >= np.abs(z),
W >= 0
All arguments must be Expression-like, and z must satisfy
z.ndim <= 1. The columns (rows) of alpha must sum to 1 when
axis=0 (axis=1).
Note: unlike PowCone3D, we make no attempt to promote
alpha to the appropriate shape. The dimensions of W and
alpha must match exactly.
Note: Dual variables are not currently implemented for this type
of constraint.
"""
_TOL_ = 1e-6
def __init__(self, W, z, alpha, axis: int = 0, constr_id=None) -> None:
Expression = cvxtypes.expression()
W = Expression.cast_to_const(W)
if not (W.is_real() and W.is_affine()):
msg = "Invalid first argument; W must be affine and real."
raise ValueError(msg)
z = Expression.cast_to_const(z)
if z.ndim > 1 or not (z.is_real() and z.is_affine()):
msg = ("Invalid second argument. z must be affine, real, "
"and have at most one z.ndim <= 1.")
raise ValueError(msg)
# Check z has one entry per cone.
if (W.ndim <= 1 and z.size > 1) or \
(W.ndim == 2 and z.size != W.shape[1-axis]) or \
(W.ndim == 1 and axis == 1):
raise ValueError(
"Argument dimensions %s and %s, with axis=%i, are incompatible."
% (W.shape, z.shape, axis))
if W.ndim == 2 and W.shape[axis] <= 1:
msg = "PowConeND requires left-hand-side to have at least two terms."
raise ValueError(msg)
alpha = Expression.cast_to_const(alpha)
if alpha.shape != W.shape:
raise ValueError("Argument dimensions %s and %s are not equal."
% (W.shape, alpha.shape))
if np.any(alpha.value <= 0):
raise ValueError("Argument alpha must be entry-wise positive.")
if np.any(np.abs(1 - np.sum(alpha.value, axis=axis)) > PowConeND._TOL_):
raise ValueError("Argument alpha must sum to 1 along axis %s." % axis)
self.W = W
self.z = z
self.alpha = alpha
self.axis = axis
if z.ndim == 0:
z = z.flatten()
super(PowConeND, self).__init__([W, z], constr_id)
def __str__(self) -> str:
return "PowND(%s, %s; %s)" % (self.W, self.z, self.alpha)
def is_imag(self) -> bool:
return False
def is_complex(self) -> bool:
return False
def get_data(self):
return [self.alpha, self.axis, self.id]
@property
def residual(self):
# TODO: The projection should be implemented directly.
from cvxpy import Minimize, Problem, Variable, hstack, norm2
if self.W.value is None or self.z.value is None:
return None
W = Variable(self.W.shape)
z = Variable(self.z.shape)
constr = [PowConeND(W, z, self.alpha, axis=self.axis)]
obj = Minimize(norm2(hstack([W.flatten(), z.flatten()]) -
hstack([self.W.flatten().value, self.z.flatten().value])))
problem = Problem(obj, constr)
return problem.solve(solver='SCS', eps=1e-8)
def num_cones(self):
return self.z.size
@property
def size(self) -> int:
cone_size = 1 + self.args[0].shape[self.axis]
return cone_size * self.num_cones()
def cone_sizes(self) -> List[int]:
cone_size = 1 + self.args[0].shape[self.axis]
return [cone_size] * self.num_cones()
[docs] def is_dcp(self, dpp: bool = False) -> bool:
"""A power cone constraint is DCP if each argument is affine.
"""
if dpp:
with scopes.dpp_scope():
args_ok = self.args[0].is_affine() and self.args[1].is_affine()
exps_ok = not isinstance(self.alpha, cvxtypes.parameter())
return args_ok and exps_ok
return True
def is_dgp(self, dpp: bool = False) -> bool:
return False
def is_dqcp(self) -> bool:
return self.is_dcp()
def save_dual_value(self, value) -> None:
dW = value[:, :-1]
dz = value[:, -1]
if self.axis == 0:
dW = dW.T
dz = dz.T
if dW.shape[1] == 1:
#NOTE: Targetting problems where duals have the shape
# (n, 1) --- dropping the extra dimension is crucial for
# the `_dual_cone` and `dual_residual` methods to work properly
dW = np.squeeze(dW)
self.dual_variables[0].save_value(dW)
self.dual_variables[1].save_value(dz)
def _dual_cone(self, *args):
"""Implements the dual cone of PowConeND See Pg 85
of the MOSEK modelling cookbook for more information"""
if args is None or args == ():
scaled_duals = self.dual_variables[0]/self.alpha
return PowConeND(scaled_duals, self.dual_variables[1], self.alpha, axis=self.axis)
else:
# some assertions for verifying `args`
args_shapes = [arg.shape for arg in args]
instance_args_shapes = [arg.shape for arg in self.args]
assert len(args) == len(self.args)
assert args_shapes == instance_args_shapes
assert args[0].value.shape == self.alpha.value.shape
return PowConeND(args[0]/self.alpha, args[1], self.alpha, axis=self.axis)