"""
UncertaintySet class: defines generic methods and attributes
of an uncertainty set in the context of robust optimization. UncertaintySet objects only
contain data which describes the set, and does not contain any Pyomo object information.
Supports the following classes of uncertainty sets:
- UncertaintySet (user defined/implemented)
- Ellipsoidal
- AxesAlignedEllipsoidal
- Polyhedral
- Box
- BudgetSet
- Cardinality/Gamma
- Discrete
- FactorModel
- IntersectedSet
"""
import abc
import functools
import math
from enum import Enum
from pyomo.common.dependencies import numpy as np, scipy as sp
from pyomo.core.base import ConcreteModel, Objective, maximize, minimize, Block
from pyomo.core.base.constraint import ConstraintList
from pyomo.core.base.var import Var, _VarData, IndexedVar
from pyomo.core.base.param import Param, _ParamData, IndexedParam
from pyomo.core.expr import value
from pyomo.opt.results import check_optimal_termination
from pyomo.contrib.pyros.util import add_bounds_for_uncertain_parameters
def uncertainty_sets(obj):
if not isinstance(obj, UncertaintySet):
raise ValueError("Expected an UncertaintySet object, instead recieved %s" % (obj,))
return obj
def column(matrix, i):
# Get column i of a given multi-dimensional list
return [row[i] for row in matrix]
class Geometry(Enum):
'''
Enum defining uncertainty set geometries
'''
LINEAR = 1
CONVEX_NONLINEAR = 2
GENERAL_NONLINEAR = 3
DISCRETE_SCENARIOS = 4
[docs]class UncertaintySet(object, metaclass=abc.ABCMeta):
'''
Base class for custom user-defined uncertainty sets.
'''
[docs] def __init__(self, **kwargs):
"""
Constructor for UncertaintySet base class
Args:
kwargs: Use the kwargs for specifying data for the UncertaintySet object. This data should be used in defining constraints in the 'set_as_constraint' function.
"""
return
@property
@abc.abstractmethod
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
raise NotImplementedError
@property
@abc.abstractmethod
def geometry(self):
"""
UncertaintySet geometry:
1 is linear,
2 is convex nonlinear,
3 is general nonlinear,
4 is discrete.
"""
raise NotImplementedError
@property
@abc.abstractmethod
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
"""
raise NotImplementedError
def is_bounded(self, config):
"""
Return True if the uncertainty set is bounded, else False.
"""
# === Determine bounds on all uncertain params
bounding_model = ConcreteModel()
bounding_model.util = Block() # So that boundedness checks work for Cardinality and FactorModel sets
bounding_model.uncertain_param_vars = IndexedVar(range(len(config.uncertain_params)), initialize=1)
for idx, param in enumerate(config.uncertain_params):
bounding_model.uncertain_param_vars[idx].set_value(
param.value, skip_validation=True)
bounding_model.add_component("uncertainty_set_constraint",
config.uncertainty_set.set_as_constraint(
uncertain_params=bounding_model.uncertain_param_vars,
model=bounding_model,
config=config
))
for idx, param in enumerate(list(bounding_model.uncertain_param_vars.values())):
bounding_model.add_component("lb_obj_" + str(idx), Objective(expr=param, sense=minimize))
bounding_model.add_component("ub_obj_" + str(idx), Objective(expr=param, sense=maximize))
for o in bounding_model.component_data_objects(Objective):
o.deactivate()
for i in range(len(bounding_model.uncertain_param_vars)):
for limit in ("lb", "ub"):
getattr(bounding_model, limit + "_obj_" + str(i)).activate()
res = config.global_solver.solve(bounding_model, tee=False)
getattr(bounding_model, limit + "_obj_" + str(i)).deactivate()
if not check_optimal_termination(res):
return False
return True
def is_nonempty(self, config):
"""
Return True if the uncertainty set is nonempty, else False.
"""
return self.is_bounded(config)
def is_valid(self, config):
"""
Return True if the uncertainty set is bounded and non-empty, else False.
"""
return self.is_nonempty(config=config) and self.is_bounded(config=config)
@abc.abstractmethod
def set_as_constraint(self, **kwargs):
"""
An uncertainty set *must* have a set_as_constraint method. UncertaintySets are instantiated with "q" as
the list of uncertain param objects. Returns a Pyomo Constraint object (could
be indexed) representing the uncertainty set for use in the separation problem
Args:
**kwargs: may be used to pass any additional information needed to generate the constraint(s)
representing the UncertaintySet
"""
pass
[docs] def point_in_set(self, point):
"""
Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False.
Args:
point: The point being checked for membership in the set.
The coordinates of the point should be supplied in the same order as the elements of ``uncertain_params``
that is to be supplied to the PyROS solve statement.
This point must match the dimension of the uncertain parameters of the set.
"""
# === Ensure point is of correct dimensionality as the uncertain parameters
if len(point) != self.dim:
raise AttributeError("Point must have same dimensions as uncertain parameters.")
m = ConcreteModel()
the_params = []
for i in range(self.dim):
m.add_component("x_%s" % i, Var(initialize=point[i]))
the_params.append(getattr(m, "x_%s" % i))
# === Generate constraint for set
set_constraint = self.set_as_constraint(uncertain_params=the_params)
# === value() returns True if the constraint is satisfied, False else.
is_in_set = all(value(con.expr) for con in set_constraint.values())
return is_in_set
@staticmethod
def add_bounds_on_uncertain_parameters(**kwargs):
"""
Numerical bounds on uncertain parameters are used in separation. This method should take a separation-type model
and update the .lb() and .ub() property for each `uncertain_param_var` member of the model to a numerical value.
This could be an inferred bound based on the uncertainty set itself, or a big-M type bound.
If the bounds need to be numerically determined, return an empty list. See PolyhedralSet and IntersectedSet as examples.
:param kwargs: the separation model and uncertainty set objects should be passed here.
:return:
"""
config = kwargs.pop('config')
model = kwargs.pop('model')
_set = config.uncertainty_set
parameter_bounds = _set.parameter_bounds
for i, p in enumerate(model.util.uncertain_param_vars.values()):
p.setlb(parameter_bounds[i][0])
p.setub(parameter_bounds[i][1])
[docs]class BoxSet(UncertaintySet):
"""
Hyper-rectangle (a.k.a. "Box")
"""
[docs] def __init__(self, bounds):
"""
BoxSet constructor
Args:
bounds: A list of tuples providing lower and upper bounds (lb, ub) for each uncertain parameter, in the same order as the 'uncertain_params' required input that is to be supplied to the PyROS solve statement.
"""
# === non-empty bounds
if len(bounds) == 0:
raise AttributeError("Vector of bounds must be non-empty")
# === Real number valued bounds
if not all(isinstance(bound, (int, float)) for tup in bounds for bound in tup):
raise AttributeError("Bounds must be real numbers.")
# === Ensure no bound is None e.g. all are bounded
if any(bound is None for tup in bounds for bound in tup):
raise AttributeError("All bounds for uncertain parameters must be real numbers, not None.")
# === Ensure each tuple has a lower and upper bound
if not all(len(b) == 2 for b in bounds):
raise AttributeError("Vector of bounds must include a finite lb and ub for each uncertain parameter")
# === Ensure each lb <= ub
if not all(bound[0] <= bound[1] for bound in bounds):
raise AttributeError("Lower bounds must be less than or equal to upper bounds")
self.bounds = bounds
self.type = "box"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return len(self.bounds)
@property
def geometry(self):
return Geometry.LINEAR
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
"""
return self.bounds
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the BoxSet uncertainty set.
Args:
uncertain_params: uncertain parameter objects for writing constraint objects
"""
conlist = ConstraintList()
conlist.construct()
set_i = list(range(len(uncertain_params)))
for i in set_i:
conlist.add(uncertain_params[i] >= self.bounds[i][0])
conlist.add(uncertain_params[i] <= self.bounds[i][1])
return conlist
[docs]class CardinalitySet(UncertaintySet):
"""
Cardinality-constrained (a.k.a "Gamma") uncertainty set
"""
[docs] def __init__(self, origin, positive_deviation, gamma):
"""
CardinalitySet constructor
Args:
origin: The origin of the set (e.g., the nominal point).
positive_deviation: Vector (``list``) of maximal deviations of each parameter.
gamma: Scalar to bound the total number of uncertain parameters that can maximally deviate from their respective 'origin'. Setting 'gamma = 0' reduces the set to the 'origin' point. Setting 'gamma' to be equal to the number of parameters produces the hyper-rectangle [origin, origin+positive_deviation]
"""
# === Real number valued data
if not all(isinstance(elem, (int, float)) for elem in origin):
raise AttributeError("Elements of origin vector must be numeric.")
if not all(isinstance(elem, (int, float)) for elem in positive_deviation):
raise AttributeError("Elements of positive_deviation vector must be numeric")
# === Dimension of positive_deviations and origin must be same
if len(origin) != len(positive_deviation):
raise AttributeError("Vectors for origin and positive_deviation must have same dimensions.")
# === Gamma between 0,1
if gamma < 0 or gamma > len(origin):
raise AttributeError("Gamma parameter must be in [0, n].")
# === positive_deviations must all be >= 0
if any(elem < 0 for elem in positive_deviation):
raise AttributeError("Elements of positive_deviations vector must be non-negative.")
# === Non-emptiness is implied
self.origin = origin
self.positive_deviation = positive_deviation
self.gamma = gamma
self.type = "cardinality"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return len(self.origin)
@property
def geometry(self):
return Geometry.LINEAR
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
"""
nom_val = self.origin
deviation = self.positive_deviation
gamma = self.gamma
parameter_bounds = [(nom_val[i], nom_val[i] + min(gamma, 1) * deviation[i]) for i in range(len(nom_val))]
return parameter_bounds
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the CardinalitySet uncertainty set.
Args:
uncertain_params: uncertain parameter objects for writing constraint objects
"""
# === Ensure dimensions
if len(uncertain_params) != len(self.origin):
raise AttributeError("Dimensions of origin and uncertain_param lists must be equal.")
model = kwargs['model']
set_i = list(range(len(uncertain_params)))
model.util.cassi = Var(set_i, initialize=0, bounds=(0, 1))
# Make n equality constraints
conlist = ConstraintList()
conlist.construct()
for i in set_i:
conlist.add(self.origin[i] + self.positive_deviation[i] * model.util.cassi[i] == uncertain_params[i])
conlist.add(sum(model.util.cassi[i] for i in set_i) <= self.gamma)
return conlist
[docs] def point_in_set(self, point):
"""
Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False.
Args:
point: the point being checked for membership in the set
"""
cassis = []
for i in range(self.dim):
if self.positive_deviation[i] > 0:
cassis.append((point[i] - self.origin[i])/self.positive_deviation[i])
if sum(cassi for cassi in cassis) <= self.gamma and \
all(cassi >= 0 and cassi <= 1 for cassi in cassis):
return True
else:
return False
[docs]class PolyhedralSet(UncertaintySet):
"""
Polyhedral uncertainty set
"""
[docs] def __init__(self, lhs_coefficients_mat, rhs_vec):
"""
PolyhedralSet constructor
Args:
lhs_coefficients_mat: Matrix of left-hand side coefficients for the linear inequality constraints defining the polyhedral set.
rhs_vec: Vector (``list``) of right-hand side values for the linear inequality constraints defining the polyhedral set.
"""
# === Real valued data
mat = np.asarray(lhs_coefficients_mat)
if not all(isinstance(elem, (int, float)) for row in lhs_coefficients_mat for elem in row):
raise AttributeError("Matrix lhs_coefficients_mat must be real-valued and numeric.")
if not all(isinstance(elem, (int, float)) for elem in rhs_vec):
raise AttributeError("Vector rhs_vec must be real-valued and numeric.")
# === Check columns of A must be same length as rhs
if mat.shape[0] != len(rhs_vec):
raise AttributeError("Rows of lhs_coefficients_mat matrix must equal length of rhs_vec list.")
# === Columns are non-zero
if mat.shape[1] == 0:
raise AttributeError("Columns of lhs_coefficients_mat must be non-zero.")
# === Matrix is not all zeros
if all(np.isclose(elem, 0) for row in lhs_coefficients_mat for elem in row):
raise AttributeError("Matrix lhs_coefficients_mat cannot be all zeroes.")
# === Non-emptiness
res = sp.optimize.linprog(c=np.zeros(mat.shape[1]), A_ub=mat, b_ub=rhs_vec, method="simplex")
if not res.success:
raise AttributeError("User-defined PolyhedralSet was determined to be empty. "
"Please check the set of constraints supplied during set construction.")
# === Boundedness
if res.status == 3:
# scipy linprog status == 3 indicates unboundedness
raise AttributeError("User-defined PolyhedralSet was determined to be unbounded. "
"Please augment the set of constraints supplied during set construction.")
self.coefficients_mat = lhs_coefficients_mat
self.rhs_vec = rhs_vec
self.type = "polyhedral"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return len(self.coefficients_mat[0])
@property
def geometry(self):
return Geometry.LINEAR
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
PolyhedralSet bounds are not computed at set construction because they cannot be algebraically determined
and require access to an optimization solver.
"""
# For the PolyhedralSet, these are numerically determined
# in the algorithm therefore they cannot presently be determined at construction of the set.
return []
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the PolyhedralSet uncertainty set.
Args:
uncertain_params: uncertain parameter objects for writing constraint objects
"""
# === Ensure valid dimensions of lhs and rhs w.r.t uncertain_params
if np.asarray(self.coefficients_mat).shape[1] != len(uncertain_params):
raise AttributeError("Columns of coefficients_mat matrix "
"must equal length of uncertain parameters list.")
set_i = list(range(len(self.coefficients_mat)))
conlist = ConstraintList()
conlist.construct()
for i in set_i:
constraint = 0
for j in range(len(uncertain_params)):
constraint += float(self.coefficients_mat[i][j]) * uncertain_params[j]
conlist.add(constraint <= float(self.rhs_vec[i]))
return conlist
@staticmethod
def add_bounds_on_uncertain_parameters(model, config):
'''
Add bounds on uncertain parameters
Args:
model: The model to add bounds on for the uncertain parameter variable objects
config: the config object for the PyROS solver instance
'''
add_bounds_for_uncertain_parameters(model=model, config=config)
return
[docs]class BudgetSet(PolyhedralSet):
"""
Budget uncertainty set
"""
[docs] def __init__(self, budget_membership_mat, rhs_vec):
"""
BudgetSet constructor
Args:
budget_membership_mat: A matrix with 0-1 entries to designate which uncertain parameters participate in each budget constraint. Here, each row is associated with a separate budget constraint.
rhs_vec: Vector (``list``) of right-hand side values for the budget constraints.
"""
# === Non-zero number of columns
mat = np.asarray(budget_membership_mat)
rhs = np.asarray(rhs_vec)
if len(mat.shape) == 1:
cols = mat.shape
else:
cols = mat.shape[1]
if cols == 0:
raise AttributeError("Budget membership matrix must have non-zero number of columns.")
# === Assert is valid matrix (same number of columns across all rows
if not all(len(row) == cols for row in budget_membership_mat):
raise AttributeError("Budget membership matrix must be a valid matrix, "
"e.g. same number of column entries across rows.")
# === Matrix dimension compatibility
if mat.shape[0] != rhs.shape[0] :
raise AttributeError("Rows of lhs_coefficients_mat matrix must equal rows of rhs_vec lists.")
# === Ensure a 0-1 matrix
if any(not np.isclose(elem, 0) and not np.isclose(elem, 1) for row in budget_membership_mat for elem in row):
raise AttributeError("Budget membership matrix must be a matrix of 0's and 1's.")
# === No all zero rows
if all(elem == 0 for row in budget_membership_mat for elem in row):
raise AttributeError("All zero rows are not permitted in the budget membership matrix.")
# === Ensure 0 <= rhs_i for all i
if any(rhs_vec[i] < 0 for i in range(len(rhs_vec))):
raise AttributeError("RHS vector entries must be >= 0.")
# === Non-emptiness is implied by the set
# === Add constraints such that uncertain params are >= 0
# === This adds >=0 bound on all parameters in the set
cols = mat.shape[1]
identity = np.identity(cols) * -1
for row in identity:
budget_membership_mat.append(row.tolist())
for i in range(identity.shape[0]):
rhs_vec.append(0)
self.coefficients_mat = budget_membership_mat
self.rhs_vec = rhs_vec
self.type = "budget"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return np.asarray(self.coefficients_mat).shape[1]
@property
def geometry(self):
return Geometry.LINEAR
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
"""
membership_mat = np.asarray(self.coefficients_mat)
rhs_vec = self.rhs_vec
parameter_bounds = []
for i in range(membership_mat.shape[1]):
col = column(membership_mat, i)
ub = min(list(col[j] * rhs_vec[j] for j in range(len(rhs_vec))))
lb = 0
parameter_bounds.append((lb, ub))
return parameter_bounds
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the BudgetSet uncertainty set.
Args:
uncertain_params: uncertain parameter objects for writing constraint objects
"""
# === Ensure matrix cols == len uncertain params
if np.asarray(self.coefficients_mat).shape[1] != len(uncertain_params):
raise AttributeError("Budget membership matrix must have compatible "
"dimensions with uncertain parameters vector.")
conlist = PolyhedralSet.set_as_constraint(self, uncertain_params)
return conlist
@staticmethod
def add_bounds_on_uncertain_parameters(model, config):
# In this case, we use the UncertaintySet class method because we have numerical parameter_bounds
UncertaintySet.add_bounds_on_uncertain_parameters(model=model, config=config)
[docs]class FactorModelSet(UncertaintySet):
"""
Factor model (a.k.a. "net-alpha" model) uncertainty set
"""
[docs] def __init__(self, origin, number_of_factors, psi_mat, beta):
"""
FactorModelSet constructor
Args:
origin: Vector (``list``) of uncertain parameter values around which deviations are restrained.
number_of_factors: Natural number representing the dimensionality of the space to which the set projects.
psi: Matrix with non-negative entries designating each uncertain parameter's contribution to each factor. Here, each row is associated with a separate uncertain parameter and each column with a separate factor.
beta: Number in [0,1] representing the fraction of the independent factors that can simultaneously attain their extreme values. Setting 'beta = 0' will enforce that as many factors will be above 0 as there will be below 0 (i.e., "zero-net-alpha" model). Setting 'beta = 1' produces the hyper-rectangle [origin - psi e, origin + psi e], where 'e' is the vector of ones.
"""
mat = np.asarray(psi_mat)
# === Numeric valued arrays
if not all(isinstance(elem, (int, float)) for elem in origin):
raise AttributeError("All elements of origin vector must be numeric.")
if not all(isinstance(elem, (int, float)) for row in psi_mat for elem in row):
raise AttributeError("All elements of psi_mat vector must be numeric.")
if not isinstance(beta, (int, float)):
raise AttributeError("Beta parameter must be numeric.")
if not isinstance(number_of_factors, (int)):
raise AttributeError("number_of_factors must be integer.")
# === Ensure dimensions of psi are n x F
if mat.shape != (len(origin), number_of_factors):
raise AttributeError("Psi matrix must be of dimensions n x F where n is dim(uncertain_params)"
"and F is number_of_factors.")
# === Ensure beta in [0,1]
if beta > 1 or beta < 0:
raise AttributeError("Beta parameter must be in [0,1].")
# === No all zero columns of psi_mat
for idx in range(mat.shape[1]):
if all(np.isclose(elem, 0) for elem in mat[:,idx]):
raise AttributeError("Psi matrix cannot have all zero columns.")
# === Psi must be strictly positive entries
for idx in range(mat.shape[1]):
if any(elem < 0 for elem in mat[:,idx]):
raise AttributeError("Psi matrix cannot have any negative entries. All factors must be non-negative.")
self.origin = origin
self.number_of_factors = number_of_factors
self.psi_mat = psi_mat
self.beta = beta
self.type = "factor_model"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return len(self.origin)
@property
def geometry(self):
return Geometry.LINEAR
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
"""
nom_val = self.origin
psi_mat = self.psi_mat
F = self.number_of_factors
beta_F = self.beta * F
floor_beta_F = math.floor(beta_F)
parameter_bounds = []
for i in range(len(nom_val)):
non_decreasing_factor_row = sorted(psi_mat[i], reverse=True)
# deviation = sum_j=1^floor(beta F) {psi_if_j} + (beta F - floor(beta F)) psi_{if_{betaF +1}}
# because indexing starts at 0, we adjust the limit on the sum and the final factor contribution
if beta_F - floor_beta_F == 0:
deviation = sum(non_decreasing_factor_row[j] for j in range(floor_beta_F - 1))
else:
deviation = sum(non_decreasing_factor_row[j] for j in range(floor_beta_F - 1)) + (
beta_F - floor_beta_F) * psi_mat[i][floor_beta_F]
lb = nom_val[i] - deviation
ub = nom_val[i] + deviation
if lb > ub:
raise AttributeError("The computed lower bound on uncertain parameters must be less than or equal to the upper bound.")
parameter_bounds.append((lb, ub))
return parameter_bounds
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the FactorModelSet uncertainty set.
Args:
uncertain_params: uncertain parameter objects for writing constraint objects
"""
model = kwargs['model']
# === Ensure dimensions
if len(uncertain_params) != len(self.origin):
raise AttributeError("Dimensions of origin and uncertain_param lists must be equal.")
# Make F-dim cassi variable
n = list(range(self.number_of_factors))
model.util.cassi = Var(n, initialize=0, bounds=(-1, 1))
conlist = ConstraintList()
conlist.construct()
disturbances = [sum(self.psi_mat[i][j] * model.util.cassi[j] for j in n)
for i in range(len(uncertain_params))]
# Make n equality constraints
for i in range(len(uncertain_params)):
conlist.add(self.origin[i] + disturbances[i] == uncertain_params[i])
conlist.add(sum(model.util.cassi[i] for i in n) <= +self.beta * self.number_of_factors)
conlist.add(sum(model.util.cassi[i] for i in n) >= -self.beta * self.number_of_factors)
return conlist
[docs] def point_in_set(self, point):
"""
Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False.
Args:
point: the point being checked for membership in the set
"""
inv_psi = np.linalg.pinv(self.psi_mat)
diff = np.asarray(list(point[i] - self.origin[i] for i in range(len(point))))
cassis = np.dot(inv_psi, np.transpose(diff))
if abs(sum(cassi for cassi in cassis)) <= self.beta * self.number_of_factors and \
all(cassi >= -1 and cassi <= 1 for cassi in cassis):
return True
else:
return False
[docs]class AxisAlignedEllipsoidalSet(UncertaintySet):
'''
Axis-aligned ellipsoidal uncertainty set
'''
[docs] def __init__(self, center, half_lengths):
"""
AxisAlignedEllipsoidalSet constructor
Args:
center: Vector (``list``) of uncertain parameter values around which deviations are restrained.
half_lengths: Vector (``list``) of half-length values representing the maximal deviations for each uncertain parameter.
"""
# === Valid data in lists
if not all(isinstance(elem, (int, float)) for elem in half_lengths):
raise AttributeError("Vector of half-lengths must be real-valued and numeric.")
if not all(isinstance(elem, (int, float)) for elem in center):
raise AttributeError("Vector center must be real-valued and numeric.")
if any(elem <= 0 for elem in half_lengths):
raise AttributeError("Half length values must be > 0.")
# === Valid variance dimensions
if not len(center) == len(half_lengths):
raise AttributeError("Half lengths and center of ellipsoid must have same dimensions.")
self.center=center
self.half_lengths=half_lengths
self.type="ellipsoidal"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return len(self.center)
@property
def geometry(self):
return Geometry.CONVEX_NONLINEAR
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
"""
nom_value = self.center
half_length =self.half_lengths
parameter_bounds = [(nom_value[i] - half_length[i], nom_value[i] + half_length[i]) for i in range(len(nom_value))]
return parameter_bounds
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the AxisAlignedEllipsoid uncertainty set.
Args:
uncertain_params: uncertain parameter objects for writing constraint objects
"""
if len(uncertain_params) != len(self.center):
raise AttributeError("Center of ellipsoid must be same dimensions as vector of uncertain parameters.")
# square and invert half lengths
inverse_squared_half_lengths = list(1.0/(a**2) for a in self.half_lengths)
# Calculate row vector of differences
diff_squared = []
# === Assume VarList uncertain_param_vars
for idx, i in enumerate(uncertain_params):
if uncertain_params[idx].is_indexed():
for index in uncertain_params[idx]:
diff_squared.append((uncertain_params[idx][index] - self.center[idx])**2)
else:
diff_squared.append((uncertain_params[idx] - self.center[idx])**2)
# Calculate inner product of difference vector and variance matrix
constraint = sum([x * y for x, y in zip(inverse_squared_half_lengths, diff_squared)])
conlist = ConstraintList()
conlist.construct()
conlist.add(constraint <= 1)
return conlist
[docs]class EllipsoidalSet(UncertaintySet):
"""
Ellipsoidal uncertainty set
"""
[docs] def __init__(self, center, shape_matrix, scale=1):
"""
EllipsoidalSet constructor
Args:
center: Vector (``list``) of uncertain parameter values around which deviations are restrained.
shape_matrix: Positive semi-definite matrix, effectively a covariance matrix for
constraint and bounds determination.
scale: Right-hand side value for the ellipsoid.
"""
# === Valid data in lists/matrixes
if not all(isinstance(elem, (int, float)) for row in shape_matrix for elem in row):
raise AttributeError("Matrix shape_matrix must be real-valued and numeric.")
if not all(isinstance(elem, (int, float)) for elem in center):
raise AttributeError("Vector center must be real-valued and numeric.")
if not isinstance(scale, (int, float)):
raise AttributeError("Ellipse scale must be a real-valued numeric.")
# === Valid matrix dimensions
num_cols = len(shape_matrix[0])
if not all(len(row) == num_cols for row in shape_matrix):
raise AttributeError("Shape matrix must have valid matrix dimensions.")
# === Ensure shape_matrix is a square matrix
array_shape_mat = np.asarray(shape_matrix)
if array_shape_mat.shape[0] != array_shape_mat.shape[1]:
raise AttributeError("Shape matrix must be square.")
# === Ensure dimensions of shape_matrix are same as dimensions of uncertain_params
if array_shape_mat.shape[1] != len(center):
raise AttributeError("Shape matrix must be "
"same dimensions as vector of uncertain parameters.")
# === Symmetric shape_matrix
if not np.all(np.abs(array_shape_mat-array_shape_mat.T) < 1e-8):
raise AttributeError("Shape matrix must be symmetric.")
# === Ensure scale is non-negative
if scale < 0:
raise AttributeError("Scale of ellipse (rhs) must be non-negative.")
# === Check if shape matrix is invertible
try:
np.linalg.inv(shape_matrix)
except np.linalg.LinAlgError as err:
raise("Error with shape matrix supplied to EllipsoidalSet object being singular. %s" % err)
# === Check is shape matrix is positive semidefinite
if not all(np.linalg.eigvals(shape_matrix) >= 0):
raise("Non positive-semidefinite shape matrix.")
# === Ensure matrix is not degenerate, for determining inferred bounds
try:
for i in range(len(shape_matrix)):
np.power(shape_matrix[i][i], 0.5)
except:
raise AttributeError("Shape matrix must be non-degenerate.")
self.center = center
self.shape_matrix = shape_matrix
self.scale = scale
self.type = "ellipsoidal"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return len(self.center)
@property
def geometry(self):
return Geometry.CONVEX_NONLINEAR
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
"""
scale = self.scale
nom_value = self.center
P = self.shape_matrix
parameter_bounds = [(nom_value[i] - np.power(P[i][i], 0.5) * scale,
nom_value[i] + np.power(P[i][i], 0.5) * scale) for i in range(self.dim)]
return parameter_bounds
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the EllipsoidalSet uncertainty set.
Args:
uncertain_params: uncertain parameter objects for writing constraint objects
"""
inv_covar = np.linalg.inv(self.shape_matrix)
if len(uncertain_params) != len(self.center):
raise AttributeError("Center of ellipsoid must be same dimensions as vector of uncertain parameters.")
# Calculate row vector of differences
diff = []
# === Assume VarList uncertain_param_vars
for idx, i in enumerate(uncertain_params):
if uncertain_params[idx].is_indexed():
for index in uncertain_params[idx]:
diff.append(uncertain_params[idx][index] - self.center[idx])
else:
diff.append(uncertain_params[idx] - self.center[idx])
# Calculate inner product of difference vector and covar matrix
product1 = [sum([x * y for x, y in zip(diff, column(inv_covar, i))]) for i in range(len(inv_covar))]
constraint = sum([x * y for x, y in zip(product1, diff)])
conlist = ConstraintList()
conlist.construct()
conlist.add(constraint <= self.scale)
return conlist
[docs]class DiscreteScenarioSet(UncertaintySet):
"""
Set of discrete scenarios (i.e., finite collection of realizations)
"""
[docs] def __init__(self, scenarios):
"""
DiscreteScenarioSet constructor
Args:
scenarios: Vector (``list``) of discrete scenarios where each scenario represents a realization of the uncertain parameters.
"""
# === Non-empty
if len(scenarios) == 0:
raise AttributeError("Scenarios list must be non-empty.")
# === Each scenario must be of real numbers
if not all(isinstance(elem, (int, float)) for d in scenarios for elem in d):
raise AttributeError("Each scenario must consist of real-number values for each parameter.")
# === Confirm all scenarios are of same dimensionality
dim = len(scenarios[0])
if not all(len(d)==dim for d in scenarios):
raise AttributeError("All points in list of scenarios must be same dimension.")
# Standardize to list of tuples
self.scenarios = list(tuple(s) for s in scenarios) # set of discrete points which are distinct realizations of uncertain params
self.type = "discrete"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return len(self.scenarios[0])
@property
def geometry(self):
return Geometry.DISCRETE_SCENARIOS
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
"""
parameter_bounds = [(min(s[i] for s in self.scenarios),
max(s[i] for s in self.scenarios)) for i in range(self.dim)]
return parameter_bounds
def is_bounded(self, config):
'''
DiscreteScenarios is bounded by default due to finiteness of the set.
:param config:
:return: True
'''
return True
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the EllipsoidalSet uncertainty set.
Args:
uncertain_params: uncertain parameter objects for writing constraint objects
"""
# === Ensure point is of correct dimensionality as the uncertain parameters
dim = len(uncertain_params)
if any(len(d) != dim for d in self.scenarios):
raise AttributeError("All scenarios must have same dimensions as uncertain parameters.")
conlist = ConstraintList()
conlist.construct()
for n in list(range(len(self.scenarios))):
for i in list(range(len(uncertain_params))):
conlist.add(uncertain_params[i] == self.scenarios[n][i])
conlist.deactivate()
return conlist
[docs] def point_in_set(self, point):
"""
Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False.
Args:
point: the point being checked for membership in the set
"""
# Round all double precision to a tolerance
num_decimals = 8
rounded_scenarios = list(list(round(num, num_decimals) for num in d) for d in self.scenarios)
rounded_point = list(round(num, num_decimals) for num in point)
return any(rounded_point==rounded_d for rounded_d in rounded_scenarios)
[docs]class IntersectionSet(UncertaintySet):
"""
Set stemming from intersecting previously constructed sets of any type
"""
[docs] def __init__(self, **kwargs):
"""
IntersectionSet constructor
Args:
**kwargs: Keyword arguments for specifying all PyROS UncertaintySet objects to be intersected.
"""
if not all(isinstance(a_set, UncertaintySet) for a_set in kwargs.values()):
raise ValueError("SetIntersection objects can only be constructed via UncertaintySet objects.")
# === dim must be defined on all UncertaintySet objects
all_sets = list(a_set for a_set in kwargs.values())
if len(all_sets) < 2:
raise AttributeError("SetIntersection requires 2 or more UncertaintySet objects.")
a_dim = all_sets[0].dim
if not all(uncertainty_set.dim == a_dim for uncertainty_set in all_sets):
raise AttributeError("Uncertainty sets being intersected must have equal dimension.")
self.all_sets = all_sets
self.type = "intersection"
@property
def dim(self):
"""
Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list.
"""
return self.all_sets[0].dim
@property
def geometry(self):
return max(self.all_sets[i].geometry.value for i in range(len(self.all_sets)))
@property
def parameter_bounds(self):
"""
Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set.
IntersectedSet bounds are not computed at set construction because they cannot be algebraically determined
and require access to an optimization solver.
"""
# For the IntersectedSet, these are numerically determined
# in the algorithm therefore they cannot presently be determined at construction of the set.
return []
[docs] def point_in_set(self, point):
"""
Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False.
Args:
point: the point being checked for membership in the set
"""
if all(a_set.point_in_set(point=point) for a_set in self.all_sets):
return True
else:
return False
def is_empty_intersection(self, uncertain_params, nlp_solver):
"""
Determine if intersection is empty
Args:
uncertain_params: list of uncertain parameters
nlp_solver: a Pyomo Solver object for solving NLPs
"""
# === Non-emptiness check for the set intersection
is_empty_intersection = True
if any(a_set.type == "discrete" for a_set in self.all_sets):
disc_sets = (a_set for a_set in self.all_sets if a_set.type == "discrete")
disc_set = min(disc_sets, key=lambda x: len(x.scenarios)) # minimum set of scenarios
# === Ensure there is at least one scenario from this discrete set which is a member of all other sets
for scenario in disc_set.scenarios:
if all(a_set.point_in_set(point=scenario) for a_set in self.all_sets):
is_empty_intersection = False
break
else:
# === Compile constraints and solve NLP
m = ConcreteModel()
m.obj = Objective(expr=0) # dummy objective required if using baron
m.param_vars = Var(uncertain_params.index_set())
for a_set in self.all_sets:
m.add_component(a_set.type + "_constraints", a_set.set_as_constraint(uncertain_params=m.param_vars))
try:
res = nlp_solver.solve(m)
except:
raise ValueError("Solver terminated with an error while checking set intersection non-emptiness.")
if check_optimal_termination(res):
is_empty_intersection = False
return is_empty_intersection
# === Define pairwise intersection function
@staticmethod
def intersect(Q1, Q2):
"""
Binary function intersecting two UncertaintySet objects
Args:
Q1: uncertainty set 1
Q2: uncertainty set 2
"""
constraints = ConstraintList()
constraints.construct()
for set in (Q1, Q2):
other = Q1 if set is Q2 else Q2
if set.type == "discrete":
intersected_scenarios = []
for point in set.scenarios:
if other.point_in_set(point=point):
intersected_scenarios.append(point)
return DiscreteScenarioSet(scenarios=intersected_scenarios)
# === This case is if both sets are continuous
return IntersectionSet(set1=Q1, set2=Q2)
return
def set_as_constraint(self, uncertain_params, **kwargs):
"""
Function to generate constraints for the IntersectedSet uncertainty set.
Args:
uncertain_params: list of uncertain param objects participating in the sets to be intersected
"""
try:
nlp_solver = kwargs["config"].global_solver
except:
raise AttributeError("set_as_constraint for SetIntersection requires access to an NLP solver via"
"the PyROS Solver config.")
is_empty_intersection = self.is_empty_intersection(uncertain_params=uncertain_params, nlp_solver=nlp_solver)
def _intersect(Q1, Q2):
return self.intersect(Q1, Q2)
if not is_empty_intersection:
Qint = functools.reduce(_intersect, self.all_sets)
if Qint.type == "discrete":
return Qint.set_as_constraint(uncertain_params=uncertain_params)
else:
conlist = ConstraintList()
conlist.construct()
for set in Qint.all_sets:
for con in list(set.set_as_constraint(uncertain_params=uncertain_params).values()):
conlist.add(con.expr)
return conlist
else:
raise AttributeError("Set intersection is empty, cannot proceed with PyROS robust optimization.")
@staticmethod
def add_bounds_on_uncertain_parameters(model, config):
"""
Add bounds on uncertain parameters
Args:
model: The model to add bounds on for the uncertain parameter variable objects
"""
add_bounds_for_uncertain_parameters(model=model, config=config)
return