Source code for psdr.domains.linineq

from __future__ import division

import numpy as np

import cvxpy as cp

from .domain import TOL
from .linquad import LinQuadDomain

from ..misc import merge

[docs]class LinIneqDomain(LinQuadDomain): r"""A domain specified by a combination of linear equality and inequality constraints. Here we build a domain specified by three kinds of constraints: bound constraints :math:`\text{lb} \le \mathbf{x} \le \text{ub}`, inequality constraints :math:`\mathbf{A} \mathbf{x} \le \mathbf{b}`, and equality constraints :math:`\mathbf{A}_{\text{eq}} \mathbf{x} = \mathbf{b}_{\text{eq}}`: .. math:: \mathcal{D} := \left \lbrace \mathbf{x} : \text{lb} \le \mathbf{x} \le \text{ub}, \ \mathbf{A} \mathbf{x} \le \mathbf{b}, \ \mathbf{A}_{\text{eq}} \mathbf{x} = \mathbf{b}_{\text{eq}} \right\rbrace \subset \mathbb{R}^m Parameters ---------- A: array-like (m,n) Matrix in left-hand side of inequality constraint b: array-like (m,) Vector in right-hand side of the ineqaluty constraint A_eq: array-like (p,n) Matrix in left-hand side of equality constraint b_eq: array-like (p,) Vector in right-hand side of equality constraint lb: array-like (n,) Vector of lower bounds ub: array-like (n,) Vector of upper bounds kwargs: dict, optional Additional parameters to pass to solvers """ def __init__(self, A = None, b = None, lb = None, ub = None, A_eq = None, b_eq = None, names = None, **kwargs): LinQuadDomain.__init__(self, A = A, b = b, lb = lb, ub = ub, A_eq = A_eq, b_eq = b_eq, names = names, **kwargs) def _isinside(self, X, tol = TOL): return self._isinside_bounds(X, tol = tol) & self._isinside_ineq(X, tol = tol) & self._isinside_eq(X, tol = tol) def _extent(self, x, p): return min(self._extent_bounds(x, p), self._extent_ineq(x, p)) def _normalized_domain(self, **kwargs): names_norm = [name + ' (normalized)' for name in self.names] return LinIneqDomain(lb = self.lb_norm, ub = self.ub_norm, A = self.A_norm, b = self.b_norm, A_eq = self.A_eq_norm, b_eq = self.b_eq_norm, names = names_norm, **merge(self.kwargs, kwargs))
[docs] def chebyshev_center(self): r"""Estimates the Chebyshev center using the constrainted least squares approach Solves the linear program finding the radius :math:`r` and Chebyshev center :math:`\mathbf{x}`. .. math:: \max_{r\in \mathbb{R}^+, \mathbf{x} \in \mathcal{D}} &\ r \\ \text{such that} & \ \mathbf{a}_i^\top \mathbf{x} + r \|\mathbf{a}_i\|_2 \le b_i where we have expressed the domain in terms of the linear inequality constraints :math:`\mathcal{D}=\lbrace \mathbf{x} : \mathbf{A}\mathbf{x} \le \mathbf{b}\rbrace` and :math:`\mathbf{a}_i^\top` are the rows of :math:`\mathbf{A}` as described in [BVNotes]_. Returns ------- center: np.ndarray(m,) Center of the domain radius: float radius of largest circle inside the domain References ---------- .. [BVNotes] https://see.stanford.edu/materials/lsocoee364a/04ConvexOptimizationProblems.pdf, page 4-19. """ m, n = self.A.shape # Merge the bound constraints into A A = [self.A] b = [self.b] for i in range(n): ei = np.zeros(n) ei[i] = 1 # Upper bound if np.isfinite(self.ub[i]): A.append(ei) b.append(self.ub[i]) # Lower bound if np.isfinite(self.lb[i]): A.append(-ei) b.append(-self.lb[i]) A = np.vstack(A) b = np.hstack(b) # See p.4-19 https://see.stanford.edu/materials/lsocoee364a/04ConvexOptimizationProblems.pdf # normA = np.sqrt( np.sum( np.power(A, 2), axis=1 ) ).reshape((A.shape[0], )) r = cp.Variable(1) x = cp.Variable(len(self)) constraints = [x.__rmatmul__(A) + normA * r <= b] if len(self.A_eq) > 0: constraints += [x.__rmatmul__(self.A_eq) == self.b_eq] problem = cp.Problem(cp.Maximize(r), constraints) problem.solve(**self.kwargs) radius = float(r.value) center = np.array(x.value).reshape(len(self)) #AA = np.hstack(( A, normA.reshape(-1,1) )) #c = np.zeros((A.shape[1]+1,)) #c[-1] = -1.0 #A_eq = np.hstack([self.A_eq, np.zeros( (self.A_eq.shape[0],1))]) #zc = linprog(c, A_ub = AA, b_ub = b, A_eq = A_eq, b_eq = self.b_eq ) #center = zc[:-1].reshape((n,)) #radius = zc[-1] #print(center) #print(self.isinside(center)) #zc = cp.Variable(len(self)+1) #prob = cp.Problem(cp.Maximize(zc[-1]), [zc.__rmatmul__(AA) <= b]) #prob.solve() #print(zc.value[0:-1]) self._radius = radius self._cheb_center = center return center, radius
@property def radius(self): try: return self._radius except: self.chebyshev_center() return self._radius @property def center(self): try: return self._cheb_center except: self.chebyshev_center() return self._cheb_center @property def Ls(self): return [] @property def ys(self): return [] @property def rhos(self): return []