1""" 2UncertaintySet class: defines generic methods and attributes 3of an uncertainty set in the context of robust optimization. UncertaintySet objects only 4contain data which describes the set, and does not contain any Pyomo object information. 5 6Supports the following classes of uncertainty sets: 7 8- UncertaintySet (user defined/implemented) 9- Ellipsoidal 10- AxesAlignedEllipsoidal 11- Polyhedral 12- Box 13- BudgetSet 14- Cardinality/Gamma 15- Discrete 16- FactorModel 17- IntersectedSet 18""" 19 20import abc 21import functools 22import math 23from enum import Enum 24from pyomo.common.dependencies import numpy as np, scipy as sp 25from pyomo.core.base import ConcreteModel, Objective, maximize, minimize, Block 26from pyomo.core.base.constraint import ConstraintList 27from pyomo.core.base.var import Var, _VarData, IndexedVar 28from pyomo.core.base.param import Param, _ParamData, IndexedParam 29from pyomo.core.expr import value 30from pyomo.opt.results import check_optimal_termination 31from pyomo.contrib.pyros.util import add_bounds_for_uncertain_parameters 32 33def uncertainty_sets(obj): 34 if not isinstance(obj, UncertaintySet): 35 raise ValueError("Expected an UncertaintySet object, instead recieved %s" % (obj,)) 36 return obj 37 38def column(matrix, i): 39 # Get column i of a given multi-dimensional list 40 return [row[i] for row in matrix] 41 42class Geometry(Enum): 43 ''' 44 Enum defining uncertainty set geometries 45 ''' 46 LINEAR = 1 47 CONVEX_NONLINEAR = 2 48 GENERAL_NONLINEAR = 3 49 DISCRETE_SCENARIOS = 4 50 51 52class UncertaintySet(object, metaclass=abc.ABCMeta): 53 ''' 54 Base class for custom user-defined uncertainty sets. 55 ''' 56 57 def __init__(self, **kwargs): 58 """ 59 Constructor for UncertaintySet base class 60 61 Args: 62 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. 63 """ 64 return 65 66 @property 67 @abc.abstractmethod 68 def dim(self): 69 """ 70 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 71 """ 72 raise NotImplementedError 73 74 @property 75 @abc.abstractmethod 76 def geometry(self): 77 """ 78 UncertaintySet geometry: 79 1 is linear, 80 2 is convex nonlinear, 81 3 is general nonlinear, 82 4 is discrete. 83 """ 84 raise NotImplementedError 85 86 @property 87 @abc.abstractmethod 88 def parameter_bounds(self): 89 """ 90 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 91 """ 92 raise NotImplementedError 93 94 def is_bounded(self, config): 95 """ 96 Return True if the uncertainty set is bounded, else False. 97 """ 98 # === Determine bounds on all uncertain params 99 bounding_model = ConcreteModel() 100 bounding_model.util = Block() # So that boundedness checks work for Cardinality and FactorModel sets 101 bounding_model.uncertain_param_vars = IndexedVar(range(len(config.uncertain_params)), initialize=1) 102 for idx, param in enumerate(config.uncertain_params): 103 bounding_model.uncertain_param_vars[idx].value = param.value 104 105 bounding_model.add_component("uncertainty_set_constraint", 106 config.uncertainty_set.set_as_constraint( 107 uncertain_params=bounding_model.uncertain_param_vars, 108 model=bounding_model, 109 config=config 110 )) 111 112 for idx, param in enumerate(list(bounding_model.uncertain_param_vars.values())): 113 bounding_model.add_component("lb_obj_" + str(idx), Objective(expr=param, sense=minimize)) 114 bounding_model.add_component("ub_obj_" + str(idx), Objective(expr=param, sense=maximize)) 115 116 for o in bounding_model.component_data_objects(Objective): 117 o.deactivate() 118 119 for i in range(len(bounding_model.uncertain_param_vars)): 120 for limit in ("lb", "ub"): 121 getattr(bounding_model, limit + "_obj_" + str(i)).activate() 122 res = config.global_solver.solve(bounding_model, tee=False) 123 getattr(bounding_model, limit + "_obj_" + str(i)).deactivate() 124 if not check_optimal_termination(res): 125 return False 126 return True 127 128 def is_nonempty(self, config): 129 """ 130 Return True if the uncertainty set is nonempty, else False. 131 """ 132 return self.is_bounded(config) 133 134 def is_valid(self, config): 135 """ 136 Return True if the uncertainty set is bounded and non-empty, else False. 137 """ 138 return self.is_nonempty(config=config) and self.is_bounded(config=config) 139 140 @abc.abstractmethod 141 def set_as_constraint(self, **kwargs): 142 """ 143 An uncertainty set *must* have a set_as_constraint method. UncertaintySets are instantiated with "q" as 144 the list of uncertain param objects. Returns a Pyomo Constraint object (could 145 be indexed) representing the uncertainty set for use in the separation problem 146 147 Args: 148 **kwargs: may be used to pass any additional information needed to generate the constraint(s) 149 representing the UncertaintySet 150 """ 151 pass 152 153 def point_in_set(self, point): 154 """ 155 Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False. 156 157 Args: 158 point: The point being checked for membership in the set. 159 The coordinates of the point should be supplied in the same order as the elements of ``uncertain_params`` 160 that is to be supplied to the PyROS solve statement. 161 This point must match the dimension of the uncertain parameters of the set. 162 """ 163 164 # === Ensure point is of correct dimensionality as the uncertain parameters 165 if len(point) != self.dim: 166 raise AttributeError("Point must have same dimensions as uncertain parameters.") 167 168 m = ConcreteModel() 169 the_params = [] 170 for i in range(self.dim): 171 m.add_component("x_%s" % i, Var(initialize=point[i])) 172 the_params.append(getattr(m, "x_%s" % i)) 173 174 # === Generate constraint for set 175 set_constraint = self.set_as_constraint(uncertain_params=the_params) 176 177 # === value() returns True if the constraint is satisfied, False else. 178 is_in_set = all(value(con.expr) for con in set_constraint.values()) 179 180 return is_in_set 181 182 @staticmethod 183 def add_bounds_on_uncertain_parameters(**kwargs): 184 """ 185 Numerical bounds on uncertain parameters are used in separation. This method should take a separation-type model 186 and update the .lb() and .ub() property for each `uncertain_param_var` member of the model to a numerical value. 187 This could be an inferred bound based on the uncertainty set itself, or a big-M type bound. 188 If the bounds need to be numerically determined, return an empty list. See PolyhedralSet and IntersectedSet as examples. 189 :param kwargs: the separation model and uncertainty set objects should be passed here. 190 :return: 191 """ 192 config = kwargs.pop('config') 193 model = kwargs.pop('model') 194 _set = config.uncertainty_set 195 parameter_bounds = _set.parameter_bounds 196 for i, p in enumerate(model.util.uncertain_param_vars.values()): 197 p.setlb(parameter_bounds[i][0]) 198 p.setub(parameter_bounds[i][1]) 199 200 201class BoxSet(UncertaintySet): 202 """ 203 Hyper-rectangle (a.k.a. "Box") 204 """ 205 206 def __init__(self, bounds): 207 """ 208 BoxSet constructor 209 210 Args: 211 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. 212 """ 213 # === non-empty bounds 214 if len(bounds) == 0: 215 raise AttributeError("Vector of bounds must be non-empty") 216 # === Real number valued bounds 217 if not all(isinstance(bound, (int, float)) for tup in bounds for bound in tup): 218 raise AttributeError("Bounds must be real numbers.") 219 # === Ensure no bound is None e.g. all are bounded 220 if any(bound is None for tup in bounds for bound in tup): 221 raise AttributeError("All bounds for uncertain parameters must be real numbers, not None.") 222 # === Ensure each tuple has a lower and upper bound 223 if not all(len(b) == 2 for b in bounds): 224 raise AttributeError("Vector of bounds must include a finite lb and ub for each uncertain parameter") 225 # === Ensure each lb <= ub 226 if not all(bound[0] <= bound[1] for bound in bounds): 227 raise AttributeError("Lower bounds must be less than or equal to upper bounds") 228 229 self.bounds = bounds 230 self.type = "box" 231 232 @property 233 def dim(self): 234 """ 235 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 236 """ 237 return len(self.bounds) 238 239 @property 240 def geometry(self): 241 return Geometry.LINEAR 242 243 @property 244 def parameter_bounds(self): 245 """ 246 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 247 """ 248 return self.bounds 249 250 def set_as_constraint(self, uncertain_params, **kwargs): 251 """ 252 Function to generate constraints for the BoxSet uncertainty set. 253 254 Args: 255 uncertain_params: uncertain parameter objects for writing constraint objects 256 """ 257 258 conlist = ConstraintList() 259 conlist.construct() 260 261 set_i = list(range(len(uncertain_params))) 262 263 for i in set_i: 264 conlist.add(uncertain_params[i] >= self.bounds[i][0]) 265 conlist.add(uncertain_params[i] <= self.bounds[i][1]) 266 267 return conlist 268 269 270class CardinalitySet(UncertaintySet): 271 """ 272 Cardinality-constrained (a.k.a "Gamma") uncertainty set 273 """ 274 275 def __init__(self, origin, positive_deviation, gamma): 276 """ 277 CardinalitySet constructor 278 279 Args: 280 origin: The origin of the set (e.g., the nominal point). 281 positive_deviation: Vector (``list``) of maximal deviations of each parameter. 282 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] 283 """ 284 # === Real number valued data 285 if not all(isinstance(elem, (int, float)) for elem in origin): 286 raise AttributeError("Elements of origin vector must be numeric.") 287 if not all(isinstance(elem, (int, float)) for elem in positive_deviation): 288 raise AttributeError("Elements of positive_deviation vector must be numeric") 289 # === Dimension of positive_deviations and origin must be same 290 if len(origin) != len(positive_deviation): 291 raise AttributeError("Vectors for origin and positive_deviation must have same dimensions.") 292 # === Gamma between 0,1 293 if gamma < 0 or gamma > len(origin): 294 raise AttributeError("Gamma parameter must be in [0, n].") 295 # === positive_deviations must all be >= 0 296 if any(elem < 0 for elem in positive_deviation): 297 raise AttributeError("Elements of positive_deviations vector must be non-negative.") 298 # === Non-emptiness is implied 299 300 self.origin = origin 301 self.positive_deviation = positive_deviation 302 self.gamma = gamma 303 self.type = "cardinality" 304 305 @property 306 def dim(self): 307 """ 308 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 309 """ 310 return len(self.origin) 311 312 @property 313 def geometry(self): 314 return Geometry.LINEAR 315 316 @property 317 def parameter_bounds(self): 318 """ 319 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 320 """ 321 322 nom_val = self.origin 323 deviation = self.positive_deviation 324 gamma = self.gamma 325 parameter_bounds = [(nom_val[i], nom_val[i] + min(gamma, 1) * deviation[i]) for i in range(len(nom_val))] 326 return parameter_bounds 327 328 def set_as_constraint(self, uncertain_params, **kwargs): 329 """ 330 Function to generate constraints for the CardinalitySet uncertainty set. 331 332 Args: 333 uncertain_params: uncertain parameter objects for writing constraint objects 334 """ 335 336 # === Ensure dimensions 337 if len(uncertain_params) != len(self.origin): 338 raise AttributeError("Dimensions of origin and uncertain_param lists must be equal.") 339 340 model = kwargs['model'] 341 set_i = list(range(len(uncertain_params))) 342 model.util.cassi = Var(set_i, initialize=0, bounds=(0, 1)) 343 344 # Make n equality constraints 345 conlist = ConstraintList() 346 conlist.construct() 347 for i in set_i: 348 conlist.add(self.origin[i] + self.positive_deviation[i] * model.util.cassi[i] == uncertain_params[i]) 349 350 conlist.add(sum(model.util.cassi[i] for i in set_i) <= self.gamma) 351 352 return conlist 353 354 def point_in_set(self, point): 355 """ 356 Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False. 357 358 Args: 359 point: the point being checked for membership in the set 360 """ 361 362 cassis = [] 363 for i in range(self.dim): 364 if self.positive_deviation[i] > 0: 365 cassis.append((point[i] - self.origin[i])/self.positive_deviation[i]) 366 367 if sum(cassi for cassi in cassis) <= self.gamma and \ 368 all(cassi >= 0 and cassi <= 1 for cassi in cassis): 369 return True 370 else: 371 return False 372 373 374class PolyhedralSet(UncertaintySet): 375 """ 376 Polyhedral uncertainty set 377 """ 378 379 def __init__(self, lhs_coefficients_mat, rhs_vec): 380 """ 381 PolyhedralSet constructor 382 383 Args: 384 lhs_coefficients_mat: Matrix of left-hand side coefficients for the linear inequality constraints defining the polyhedral set. 385 rhs_vec: Vector (``list``) of right-hand side values for the linear inequality constraints defining the polyhedral set. 386 """ 387 388 # === Real valued data 389 mat = np.asarray(lhs_coefficients_mat) 390 if not all(isinstance(elem, (int, float)) for row in lhs_coefficients_mat for elem in row): 391 raise AttributeError("Matrix lhs_coefficients_mat must be real-valued and numeric.") 392 if not all(isinstance(elem, (int, float)) for elem in rhs_vec): 393 raise AttributeError("Vector rhs_vec must be real-valued and numeric.") 394 # === Check columns of A must be same length as rhs 395 if mat.shape[0] != len(rhs_vec): 396 raise AttributeError("Rows of lhs_coefficients_mat matrix must equal length of rhs_vec list.") 397 # === Columns are non-zero 398 if mat.shape[1] == 0: 399 raise AttributeError("Columns of lhs_coefficients_mat must be non-zero.") 400 # === Matrix is not all zeros 401 if all(np.isclose(elem, 0) for row in lhs_coefficients_mat for elem in row): 402 raise AttributeError("Matrix lhs_coefficients_mat cannot be all zeroes.") 403 # === Non-emptiness 404 res = sp.optimize.linprog(c=np.zeros(mat.shape[1]), A_ub=mat, b_ub=rhs_vec, method="simplex") 405 if not res.success: 406 raise AttributeError("User-defined PolyhedralSet was determined to be empty. " 407 "Please check the set of constraints supplied during set construction.") 408 # === Boundedness 409 if res.status == 3: 410 # scipy linprog status == 3 indicates unboundedness 411 raise AttributeError("User-defined PolyhedralSet was determined to be unbounded. " 412 "Please augment the set of constraints supplied during set construction.") 413 414 415 self.coefficients_mat = lhs_coefficients_mat 416 self.rhs_vec = rhs_vec 417 self.type = "polyhedral" 418 419 @property 420 def dim(self): 421 """ 422 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 423 """ 424 return len(self.coefficients_mat[0]) 425 426 @property 427 def geometry(self): 428 return Geometry.LINEAR 429 430 @property 431 def parameter_bounds(self): 432 """ 433 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 434 PolyhedralSet bounds are not computed at set construction because they cannot be algebraically determined 435 and require access to an optimization solver. 436 """ 437 # For the PolyhedralSet, these are numerically determined 438 # in the algorithm therefore they cannot presently be determined at construction of the set. 439 return [] 440 441 def set_as_constraint(self, uncertain_params, **kwargs): 442 """ 443 Function to generate constraints for the PolyhedralSet uncertainty set. 444 445 Args: 446 uncertain_params: uncertain parameter objects for writing constraint objects 447 """ 448 449 # === Ensure valid dimensions of lhs and rhs w.r.t uncertain_params 450 if np.asarray(self.coefficients_mat).shape[1] != len(uncertain_params): 451 raise AttributeError("Columns of coefficients_mat matrix " 452 "must equal length of uncertain parameters list.") 453 454 set_i = list(range(len(self.coefficients_mat))) 455 456 conlist = ConstraintList() 457 conlist.construct() 458 459 for i in set_i: 460 constraint = 0 461 for j in range(len(uncertain_params)): 462 constraint += float(self.coefficients_mat[i][j]) * uncertain_params[j] 463 conlist.add(constraint <= float(self.rhs_vec[i])) 464 465 return conlist 466 467 @staticmethod 468 def add_bounds_on_uncertain_parameters(model, config): 469 ''' 470 Add bounds on uncertain parameters 471 472 Args: 473 model: The model to add bounds on for the uncertain parameter variable objects 474 config: the config object for the PyROS solver instance 475 ''' 476 add_bounds_for_uncertain_parameters(model=model, config=config) 477 return 478 479 480class BudgetSet(PolyhedralSet): 481 """ 482 Budget uncertainty set 483 """ 484 485 def __init__(self, budget_membership_mat, rhs_vec): 486 """ 487 BudgetSet constructor 488 489 Args: 490 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. 491 rhs_vec: Vector (``list``) of right-hand side values for the budget constraints. 492 """ 493 # === Non-zero number of columns 494 mat = np.asarray(budget_membership_mat) 495 rhs = np.asarray(rhs_vec) 496 497 if len(mat.shape) == 1: 498 cols = mat.shape 499 else: 500 cols = mat.shape[1] 501 if cols == 0: 502 raise AttributeError("Budget membership matrix must have non-zero number of columns.") 503 # === Assert is valid matrix (same number of columns across all rows 504 if not all(len(row) == cols for row in budget_membership_mat): 505 raise AttributeError("Budget membership matrix must be a valid matrix, " 506 "e.g. same number of column entries across rows.") 507 # === Matrix dimension compatibility 508 if mat.shape[0] != rhs.shape[0] : 509 raise AttributeError("Rows of lhs_coefficients_mat matrix must equal rows of rhs_vec lists.") 510 # === Ensure a 0-1 matrix 511 if any(not np.isclose(elem, 0) and not np.isclose(elem, 1) for row in budget_membership_mat for elem in row): 512 raise AttributeError("Budget membership matrix must be a matrix of 0's and 1's.") 513 # === No all zero rows 514 if all(elem == 0 for row in budget_membership_mat for elem in row): 515 raise AttributeError("All zero rows are not permitted in the budget membership matrix.") 516 517 # === Ensure 0 <= rhs_i for all i 518 if any(rhs_vec[i] < 0 for i in range(len(rhs_vec))): 519 raise AttributeError("RHS vector entries must be >= 0.") 520 # === Non-emptiness is implied by the set 521 522 # === Add constraints such that uncertain params are >= 0 523 # === This adds >=0 bound on all parameters in the set 524 cols = mat.shape[1] 525 identity = np.identity(cols) * -1 526 for row in identity: 527 budget_membership_mat.append(row.tolist()) 528 529 for i in range(identity.shape[0]): 530 rhs_vec.append(0) 531 532 self.coefficients_mat = budget_membership_mat 533 self.rhs_vec = rhs_vec 534 535 self.type = "budget" 536 537 @property 538 def dim(self): 539 """ 540 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 541 """ 542 return np.asarray(self.coefficients_mat).shape[1] 543 544 @property 545 def geometry(self): 546 return Geometry.LINEAR 547 548 @property 549 def parameter_bounds(self): 550 """ 551 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 552 """ 553 membership_mat = np.asarray(self.coefficients_mat) 554 rhs_vec = self.rhs_vec 555 parameter_bounds = [] 556 for i in range(membership_mat.shape[1]): 557 col = column(membership_mat, i) 558 ub = min(list(col[j] * rhs_vec[j] for j in range(len(rhs_vec)))) 559 lb = 0 560 parameter_bounds.append((lb, ub)) 561 return parameter_bounds 562 563 def set_as_constraint(self, uncertain_params, **kwargs): 564 """ 565 Function to generate constraints for the BudgetSet uncertainty set. 566 567 Args: 568 uncertain_params: uncertain parameter objects for writing constraint objects 569 """ 570 571 # === Ensure matrix cols == len uncertain params 572 if np.asarray(self.coefficients_mat).shape[1] != len(uncertain_params): 573 raise AttributeError("Budget membership matrix must have compatible " 574 "dimensions with uncertain parameters vector.") 575 576 conlist = PolyhedralSet.set_as_constraint(self, uncertain_params) 577 return conlist 578 579 @staticmethod 580 def add_bounds_on_uncertain_parameters(model, config): 581 # In this case, we use the UncertaintySet class method because we have numerical parameter_bounds 582 UncertaintySet.add_bounds_on_uncertain_parameters(model=model, config=config) 583 584class FactorModelSet(UncertaintySet): 585 """ 586 Factor model (a.k.a. "net-alpha" model) uncertainty set 587 """ 588 589 def __init__(self, origin, number_of_factors, psi_mat, beta): 590 """ 591 FactorModelSet constructor 592 593 Args: 594 origin: Vector (``list``) of uncertain parameter values around which deviations are restrained. 595 number_of_factors: Natural number representing the dimensionality of the space to which the set projects. 596 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. 597 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. 598 """ 599 mat = np.asarray(psi_mat) 600 # === Numeric valued arrays 601 if not all(isinstance(elem, (int, float)) for elem in origin): 602 raise AttributeError("All elements of origin vector must be numeric.") 603 if not all(isinstance(elem, (int, float)) for row in psi_mat for elem in row): 604 raise AttributeError("All elements of psi_mat vector must be numeric.") 605 if not isinstance(beta, (int, float)): 606 raise AttributeError("Beta parameter must be numeric.") 607 if not isinstance(number_of_factors, (int)): 608 raise AttributeError("number_of_factors must be integer.") 609 # === Ensure dimensions of psi are n x F 610 if mat.shape != (len(origin), number_of_factors): 611 raise AttributeError("Psi matrix must be of dimensions n x F where n is dim(uncertain_params)" 612 "and F is number_of_factors.") 613 # === Ensure beta in [0,1] 614 if beta > 1 or beta < 0: 615 raise AttributeError("Beta parameter must be in [0,1].") 616 # === No all zero columns of psi_mat 617 for idx in range(mat.shape[1]): 618 if all(np.isclose(elem, 0) for elem in mat[:,idx]): 619 raise AttributeError("Psi matrix cannot have all zero columns.") 620 # === Psi must be strictly positive entries 621 for idx in range(mat.shape[1]): 622 if any(elem < 0 for elem in mat[:,idx]): 623 raise AttributeError("Psi matrix cannot have any negative entries. All factors must be non-negative.") 624 625 self.origin = origin 626 self.number_of_factors = number_of_factors 627 self.psi_mat = psi_mat 628 self.beta = beta 629 self.type = "factor_model" 630 631 @property 632 def dim(self): 633 """ 634 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 635 """ 636 return len(self.origin) 637 638 @property 639 def geometry(self): 640 return Geometry.LINEAR 641 642 @property 643 def parameter_bounds(self): 644 """ 645 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 646 """ 647 nom_val = self.origin 648 psi_mat = self.psi_mat 649 650 F = self.number_of_factors 651 beta_F = self.beta * F 652 floor_beta_F = math.floor(beta_F) 653 parameter_bounds = [] 654 for i in range(len(nom_val)): 655 non_decreasing_factor_row = sorted(psi_mat[i], reverse=True) 656 # deviation = sum_j=1^floor(beta F) {psi_if_j} + (beta F - floor(beta F)) psi_{if_{betaF +1}} 657 # because indexing starts at 0, we adjust the limit on the sum and the final factor contribution 658 if beta_F - floor_beta_F == 0: 659 deviation = sum(non_decreasing_factor_row[j] for j in range(floor_beta_F - 1)) 660 else: 661 deviation = sum(non_decreasing_factor_row[j] for j in range(floor_beta_F - 1)) + ( 662 beta_F - floor_beta_F) * psi_mat[i][floor_beta_F] 663 lb = nom_val[i] - deviation 664 ub = nom_val[i] + deviation 665 if lb > ub: 666 raise AttributeError("The computed lower bound on uncertain parameters must be less than or equal to the upper bound.") 667 parameter_bounds.append((lb, ub)) 668 return parameter_bounds 669 670 def set_as_constraint(self, uncertain_params, **kwargs): 671 """ 672 Function to generate constraints for the FactorModelSet uncertainty set. 673 674 Args: 675 uncertain_params: uncertain parameter objects for writing constraint objects 676 """ 677 model = kwargs['model'] 678 679 # === Ensure dimensions 680 if len(uncertain_params) != len(self.origin): 681 raise AttributeError("Dimensions of origin and uncertain_param lists must be equal.") 682 683 # Make F-dim cassi variable 684 n = list(range(self.number_of_factors)) 685 model.util.cassi = Var(n, initialize=0, bounds=(-1, 1)) 686 687 conlist = ConstraintList() 688 conlist.construct() 689 690 disturbances = [sum(self.psi_mat[i][j] * model.util.cassi[j] for j in n) 691 for i in range(len(uncertain_params))] 692 693 # Make n equality constraints 694 for i in range(len(uncertain_params)): 695 conlist.add(self.origin[i] + disturbances[i] == uncertain_params[i]) 696 conlist.add(sum(model.util.cassi[i] for i in n) <= +self.beta * self.number_of_factors) 697 conlist.add(sum(model.util.cassi[i] for i in n) >= -self.beta * self.number_of_factors) 698 return conlist 699 700 701 def point_in_set(self, point): 702 """ 703 Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False. 704 705 Args: 706 point: the point being checked for membership in the set 707 """ 708 inv_psi = np.linalg.pinv(self.psi_mat) 709 diff = np.asarray(list(point[i] - self.origin[i] for i in range(len(point)))) 710 cassis = np.dot(inv_psi, np.transpose(diff)) 711 712 if abs(sum(cassi for cassi in cassis)) <= self.beta * self.number_of_factors and \ 713 all(cassi >= -1 and cassi <= 1 for cassi in cassis): 714 return True 715 else: 716 return False 717 718 719class AxisAlignedEllipsoidalSet(UncertaintySet): 720 ''' 721 Axis-aligned ellipsoidal uncertainty set 722 ''' 723 def __init__(self, center, half_lengths): 724 """ 725 AxisAlignedEllipsoidalSet constructor 726 727 Args: 728 center: Vector (``list``) of uncertain parameter values around which deviations are restrained. 729 half_lengths: Vector (``list``) of half-length values representing the maximal deviations for each uncertain parameter. 730 """ 731 # === Valid data in lists 732 if not all(isinstance(elem, (int, float)) for elem in half_lengths): 733 raise AttributeError("Vector of half-lengths must be real-valued and numeric.") 734 if not all(isinstance(elem, (int, float)) for elem in center): 735 raise AttributeError("Vector center must be real-valued and numeric.") 736 if any(elem <= 0 for elem in half_lengths): 737 raise AttributeError("Half length values must be > 0.") 738 # === Valid variance dimensions 739 if not len(center) == len(half_lengths): 740 raise AttributeError("Half lengths and center of ellipsoid must have same dimensions.") 741 742 self.center=center 743 self.half_lengths=half_lengths 744 self.type="ellipsoidal" 745 746 @property 747 def dim(self): 748 """ 749 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 750 """ 751 return len(self.center) 752 753 @property 754 def geometry(self): 755 return Geometry.CONVEX_NONLINEAR 756 757 @property 758 def parameter_bounds(self): 759 """ 760 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 761 """ 762 nom_value = self.center 763 half_length =self.half_lengths 764 parameter_bounds = [(nom_value[i] - half_length[i], nom_value[i] + half_length[i]) for i in range(len(nom_value))] 765 return parameter_bounds 766 767 def set_as_constraint(self, uncertain_params, **kwargs): 768 """ 769 Function to generate constraints for the AxisAlignedEllipsoid uncertainty set. 770 771 Args: 772 uncertain_params: uncertain parameter objects for writing constraint objects 773 """ 774 if len(uncertain_params) != len(self.center): 775 raise AttributeError("Center of ellipsoid must be same dimensions as vector of uncertain parameters.") 776 # square and invert half lengths 777 inverse_squared_half_lengths = list(1.0/(a**2) for a in self.half_lengths) 778 # Calculate row vector of differences 779 diff_squared = [] 780 # === Assume VarList uncertain_param_vars 781 for idx, i in enumerate(uncertain_params): 782 if uncertain_params[idx].is_indexed(): 783 for index in uncertain_params[idx]: 784 diff_squared.append((uncertain_params[idx][index] - self.center[idx])**2) 785 else: 786 diff_squared.append((uncertain_params[idx] - self.center[idx])**2) 787 788 # Calculate inner product of difference vector and variance matrix 789 constraint = sum([x * y for x, y in zip(inverse_squared_half_lengths, diff_squared)]) 790 791 conlist = ConstraintList() 792 conlist.construct() 793 conlist.add(constraint <= 1) 794 return conlist 795 796 797class EllipsoidalSet(UncertaintySet): 798 """ 799 Ellipsoidal uncertainty set 800 """ 801 802 def __init__(self, center, shape_matrix, scale=1): 803 """ 804 EllipsoidalSet constructor 805 806 Args: 807 center: Vector (``list``) of uncertain parameter values around which deviations are restrained. 808 shape_matrix: Positive semi-definite matrix, effectively a covariance matrix for 809 constraint and bounds determination. 810 scale: Right-hand side value for the ellipsoid. 811 """ 812 813 # === Valid data in lists/matrixes 814 if not all(isinstance(elem, (int, float)) for row in shape_matrix for elem in row): 815 raise AttributeError("Matrix shape_matrix must be real-valued and numeric.") 816 if not all(isinstance(elem, (int, float)) for elem in center): 817 raise AttributeError("Vector center must be real-valued and numeric.") 818 if not isinstance(scale, (int, float)): 819 raise AttributeError("Ellipse scale must be a real-valued numeric.") 820 # === Valid matrix dimensions 821 num_cols = len(shape_matrix[0]) 822 if not all(len(row) == num_cols for row in shape_matrix): 823 raise AttributeError("Shape matrix must have valid matrix dimensions.") 824 # === Ensure shape_matrix is a square matrix 825 array_shape_mat = np.asarray(shape_matrix) 826 if array_shape_mat.shape[0] != array_shape_mat.shape[1]: 827 raise AttributeError("Shape matrix must be square.") 828 # === Ensure dimensions of shape_matrix are same as dimensions of uncertain_params 829 if array_shape_mat.shape[1] != len(center): 830 raise AttributeError("Shape matrix must be " 831 "same dimensions as vector of uncertain parameters.") 832 # === Symmetric shape_matrix 833 if not np.all(np.abs(array_shape_mat-array_shape_mat.T) < 1e-8): 834 raise AttributeError("Shape matrix must be symmetric.") 835 # === Ensure scale is non-negative 836 if scale < 0: 837 raise AttributeError("Scale of ellipse (rhs) must be non-negative.") 838 # === Check if shape matrix is invertible 839 try: 840 np.linalg.inv(shape_matrix) 841 except np.linalg.LinAlgError as err: 842 raise("Error with shape matrix supplied to EllipsoidalSet object being singular. %s" % err) 843 # === Check is shape matrix is positive semidefinite 844 if not all(np.linalg.eigvals(shape_matrix) >= 0): 845 raise("Non positive-semidefinite shape matrix.") 846 # === Ensure matrix is not degenerate, for determining inferred bounds 847 try: 848 for i in range(len(shape_matrix)): 849 np.power(shape_matrix[i][i], 0.5) 850 except: 851 raise AttributeError("Shape matrix must be non-degenerate.") 852 853 self.center = center 854 self.shape_matrix = shape_matrix 855 self.scale = scale 856 self.type = "ellipsoidal" 857 858 @property 859 def dim(self): 860 """ 861 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 862 """ 863 return len(self.center) 864 865 @property 866 def geometry(self): 867 return Geometry.CONVEX_NONLINEAR 868 869 @property 870 def parameter_bounds(self): 871 """ 872 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 873 """ 874 scale = self.scale 875 nom_value = self.center 876 P = self.shape_matrix 877 parameter_bounds = [(nom_value[i] - np.power(P[i][i], 0.5) * scale, 878 nom_value[i] + np.power(P[i][i], 0.5) * scale) for i in range(self.dim)] 879 return parameter_bounds 880 881 def set_as_constraint(self, uncertain_params, **kwargs): 882 """ 883 Function to generate constraints for the EllipsoidalSet uncertainty set. 884 885 Args: 886 uncertain_params: uncertain parameter objects for writing constraint objects 887 """ 888 inv_covar = np.linalg.inv(self.shape_matrix) 889 890 if len(uncertain_params) != len(self.center): 891 raise AttributeError("Center of ellipsoid must be same dimensions as vector of uncertain parameters.") 892 893 # Calculate row vector of differences 894 diff = [] 895 # === Assume VarList uncertain_param_vars 896 for idx, i in enumerate(uncertain_params): 897 if uncertain_params[idx].is_indexed(): 898 for index in uncertain_params[idx]: 899 diff.append(uncertain_params[idx][index] - self.center[idx]) 900 else: 901 diff.append(uncertain_params[idx] - self.center[idx]) 902 903 # Calculate inner product of difference vector and covar matrix 904 product1 = [sum([x * y for x, y in zip(diff, column(inv_covar, i))]) for i in range(len(inv_covar))] 905 constraint = sum([x * y for x, y in zip(product1, diff)]) 906 907 conlist = ConstraintList() 908 conlist.construct() 909 conlist.add(constraint <= self.scale) 910 return conlist 911 912 913class DiscreteScenarioSet(UncertaintySet): 914 """ 915 Set of discrete scenarios (i.e., finite collection of realizations) 916 """ 917 918 def __init__(self, scenarios): 919 """ 920 DiscreteScenarioSet constructor 921 922 Args: 923 scenarios: Vector (``list``) of discrete scenarios where each scenario represents a realization of the uncertain parameters. 924 """ 925 926 # === Non-empty 927 if len(scenarios) == 0: 928 raise AttributeError("Scenarios list must be non-empty.") 929 # === Each scenario must be of real numbers 930 if not all(isinstance(elem, (int, float)) for d in scenarios for elem in d): 931 raise AttributeError("Each scenario must consist of real-number values for each parameter.") 932 # === Confirm all scenarios are of same dimensionality 933 dim = len(scenarios[0]) 934 if not all(len(d)==dim for d in scenarios): 935 raise AttributeError("All points in list of scenarios must be same dimension.") 936 937 # Standardize to list of tuples 938 self.scenarios = list(tuple(s) for s in scenarios) # set of discrete points which are distinct realizations of uncertain params 939 self.type = "discrete" 940 941 @property 942 def dim(self): 943 """ 944 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 945 """ 946 return len(self.scenarios[0]) 947 948 @property 949 def geometry(self): 950 return Geometry.DISCRETE_SCENARIOS 951 952 @property 953 def parameter_bounds(self): 954 """ 955 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 956 """ 957 parameter_bounds = [(min(s[i] for s in self.scenarios), 958 max(s[i] for s in self.scenarios)) for i in range(self.dim)] 959 return parameter_bounds 960 961 def is_bounded(self, config): 962 ''' 963 DiscreteScenarios is bounded by default due to finiteness of the set. 964 :param config: 965 :return: True 966 ''' 967 return True 968 969 def set_as_constraint(self, uncertain_params, **kwargs): 970 """ 971 Function to generate constraints for the EllipsoidalSet uncertainty set. 972 973 Args: 974 uncertain_params: uncertain parameter objects for writing constraint objects 975 """ 976 # === Ensure point is of correct dimensionality as the uncertain parameters 977 dim = len(uncertain_params) 978 if any(len(d) != dim for d in self.scenarios): 979 raise AttributeError("All scenarios must have same dimensions as uncertain parameters.") 980 981 conlist = ConstraintList() 982 conlist.construct() 983 984 for n in list(range(len(self.scenarios))): 985 for i in list(range(len(uncertain_params))): 986 conlist.add(uncertain_params[i] == self.scenarios[n][i]) 987 988 conlist.deactivate() 989 return conlist 990 991 def point_in_set(self, point): 992 """ 993 Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False. 994 995 Args: 996 point: the point being checked for membership in the set 997 """ 998 # Round all double precision to a tolerance 999 num_decimals = 8 1000 rounded_scenarios = list(list(round(num, num_decimals) for num in d) for d in self.scenarios) 1001 rounded_point = list(round(num, num_decimals) for num in point) 1002 1003 return any(rounded_point==rounded_d for rounded_d in rounded_scenarios) 1004 1005 1006class IntersectionSet(UncertaintySet): 1007 """ 1008 Set stemming from intersecting previously constructed sets of any type 1009 """ 1010 1011 def __init__(self, **kwargs): 1012 """ 1013 IntersectionSet constructor 1014 1015 Args: 1016 **kwargs: Keyword arguments for specifying all PyROS UncertaintySet objects to be intersected. 1017 """ 1018 if not all(isinstance(a_set, UncertaintySet) for a_set in kwargs.values()): 1019 raise ValueError("SetIntersection objects can only be constructed via UncertaintySet objects.") 1020 1021 # === dim must be defined on all UncertaintySet objects 1022 all_sets = list(a_set for a_set in kwargs.values()) 1023 if len(all_sets) < 2: 1024 raise AttributeError("SetIntersection requires 2 or more UncertaintySet objects.") 1025 1026 a_dim = all_sets[0].dim 1027 if not all(uncertainty_set.dim == a_dim for uncertainty_set in all_sets): 1028 raise AttributeError("Uncertainty sets being intersected must have equal dimension.") 1029 1030 self.all_sets = all_sets 1031 self.type = "intersection" 1032 1033 @property 1034 def dim(self): 1035 """ 1036 Dimension of the uncertainty set, i.e., number of parameters in “uncertain_params” list. 1037 """ 1038 return self.all_sets[0].dim 1039 1040 @property 1041 def geometry(self): 1042 return max(self.all_sets[i].geometry.value for i in range(len(self.all_sets))) 1043 1044 @property 1045 def parameter_bounds(self): 1046 """ 1047 Bounds on the realizations of the uncertain parameters, as inferred from the uncertainty set. 1048 IntersectedSet bounds are not computed at set construction because they cannot be algebraically determined 1049 and require access to an optimization solver. 1050 """ 1051 # For the IntersectedSet, these are numerically determined 1052 # in the algorithm therefore they cannot presently be determined at construction of the set. 1053 return [] 1054 1055 def point_in_set(self, point): 1056 """ 1057 Calculates if supplied ``point`` is contained in the uncertainty set. Returns True or False. 1058 1059 Args: 1060 point: the point being checked for membership in the set 1061 """ 1062 if all(a_set.point_in_set(point=point) for a_set in self.all_sets): 1063 return True 1064 else: 1065 return False 1066 1067 def is_empty_intersection(self, uncertain_params, nlp_solver): 1068 """ 1069 Determine if intersection is empty 1070 1071 Args: 1072 uncertain_params: list of uncertain parameters 1073 nlp_solver: a Pyomo Solver object for solving NLPs 1074 """ 1075 1076 # === Non-emptiness check for the set intersection 1077 is_empty_intersection = True 1078 if any(a_set.type == "discrete" for a_set in self.all_sets): 1079 disc_sets = (a_set for a_set in self.all_sets if a_set.type == "discrete") 1080 disc_set = min(disc_sets, key=lambda x: len(x.scenarios)) # minimum set of scenarios 1081 # === Ensure there is at least one scenario from this discrete set which is a member of all other sets 1082 for scenario in disc_set.scenarios: 1083 if all(a_set.point_in_set(point=scenario) for a_set in self.all_sets): 1084 is_empty_intersection = False 1085 break 1086 else: 1087 # === Compile constraints and solve NLP 1088 m = ConcreteModel() 1089 m.obj = Objective(expr=0) # dummy objective required if using baron 1090 m.param_vars = Var(uncertain_params.index_set()) 1091 for a_set in self.all_sets: 1092 m.add_component(a_set.type + "_constraints", a_set.set_as_constraint(uncertain_params=m.param_vars)) 1093 try: 1094 res = nlp_solver.solve(m) 1095 except: 1096 raise ValueError("Solver terminated with an error while checking set intersection non-emptiness.") 1097 if check_optimal_termination(res): 1098 is_empty_intersection = False 1099 return is_empty_intersection 1100 1101 # === Define pairwise intersection function 1102 @staticmethod 1103 def intersect(Q1, Q2): 1104 """ 1105 Binary function intersecting two UncertaintySet objects 1106 Args: 1107 Q1: uncertainty set 1 1108 Q2: uncertainty set 2 1109 """ 1110 constraints = ConstraintList() 1111 constraints.construct() 1112 1113 for set in (Q1, Q2): 1114 other = Q1 if set is Q2 else Q2 1115 if set.type == "discrete": 1116 intersected_scenarios = [] 1117 for point in set.scenarios: 1118 if other.point_in_set(point=point): 1119 intersected_scenarios.append(point) 1120 return DiscreteScenarioSet(scenarios=intersected_scenarios) 1121 1122 # === This case is if both sets are continuous 1123 return IntersectionSet(set1=Q1, set2=Q2) 1124 1125 return 1126 1127 def set_as_constraint(self, uncertain_params, **kwargs): 1128 """ 1129 Function to generate constraints for the IntersectedSet uncertainty set. 1130 Args: 1131 uncertain_params: list of uncertain param objects participating in the sets to be intersected 1132 """ 1133 try: 1134 nlp_solver = kwargs["config"].global_solver 1135 except: 1136 raise AttributeError("set_as_constraint for SetIntersection requires access to an NLP solver via" 1137 "the PyROS Solver config.") 1138 is_empty_intersection = self.is_empty_intersection(uncertain_params=uncertain_params, nlp_solver=nlp_solver) 1139 1140 def _intersect(Q1, Q2): 1141 return self.intersect(Q1, Q2) 1142 1143 if not is_empty_intersection: 1144 Qint = functools.reduce(_intersect, self.all_sets) 1145 1146 if Qint.type == "discrete": 1147 return Qint.set_as_constraint(uncertain_params=uncertain_params) 1148 else: 1149 conlist = ConstraintList() 1150 conlist.construct() 1151 for set in Qint.all_sets: 1152 for con in list(set.set_as_constraint(uncertain_params=uncertain_params).values()): 1153 conlist.add(con.expr) 1154 return conlist 1155 else: 1156 raise AttributeError("Set intersection is empty, cannot proceed with PyROS robust optimization.") 1157 1158 @staticmethod 1159 def add_bounds_on_uncertain_parameters(model, config): 1160 """ 1161 Add bounds on uncertain parameters 1162 1163 Args: 1164 model: The model to add bounds on for the uncertain parameter variable objects 1165 """ 1166 1167 add_bounds_for_uncertain_parameters(model=model, config=config) 1168 return 1169