diff --git a/romodel/reformulate.py b/romodel/reformulate.py deleted file mode 100644 index 4cf9b53..0000000 --- a/romodel/reformulate.py +++ /dev/null @@ -1,564 +0,0 @@ -from pyomo.environ import (Constraint, - Var, - quicksum, - sqrt, - Objective, - maximize, - minimize, - ConstraintList, - NonNegativeReals, - NonPositiveReals, - Block, - native_numeric_types) -from pyomo.environ import sqrt as pyomo_sqrt -from pyomo.core import Transformation, TransformationFactory -from pyomo.repn import generate_standard_repn -from romodel.visitor import _expression_is_uncertain -from romodel.generator import RobustConstraint -from romodel.duality import create_linear_dual_from -from romodel.uncset import WarpedGPSet, GPSet -from itertools import chain -from romodel.util import collect_uncparam -from pyomo.core.expr.visitor import replace_expressions -import numpy as np - - -class BaseRobustTransformation(Transformation): - def __init__(self): - self._fixed_unc_params = [] - self._fixed_components = {} - - def fix_component(self, instance, component=Var): - fixed = [] - for c in instance.component_data_objects(component, active=True): - if not c.is_fixed(): - c.fix() - fixed.append(c) - self._fixed_components[component] = fixed - - def unfix_component(self, component=Var): - for c in self._fixed_components[component]: - c.unfix() - self._fixed_components[component] = None - - def get_uncertain_components(self, instance, component=Constraint): - """ Return all uncertain components of type `component`. """ - comp_list = instance.component_data_objects(component, active=True) - for c in comp_list: - if hasattr(c, 'body'): - # Constraints - expr = c.body - else: - # Objective - expr = c.expr - - if _expression_is_uncertain(expr): - yield c - - def generate_repn_param(self, instance, cdata): - self.fix_component(instance, component=Var) - if hasattr(cdata, 'body'): - expr = cdata.body - else: - expr = cdata.expr - repn = generate_standard_repn(expr, compute_values=False) - self.unfix_component(component=Var) - return repn - - -@TransformationFactory.register('romodel.ellipsoidal', - doc="Ellipsoidal Counterpart") -class EllipsoidalTransformation(BaseRobustTransformation): - def _apply_to(self, instance, root=False): - for c in chain(self.get_uncertain_components(instance), - self.get_uncertain_components(instance, - component=Objective)): - # Collect unc. parameter and unc. set - param = collect_uncparam(c) - uncset = param._uncset - assert uncset is not None, ("No uncertainty set provided for " - "uncertain parameter {}." - .format(param.name)) - - # Check if uncertainty set is empty - assert not uncset.is_empty(), ("{} does not have any " - "constraints.".format(uncset.name)) - # Check if uncertainty set is ellipsoidal - if not uncset.is_ellipsoidal(): - continue - - if c.ctype is Constraint: - assert not c.equality, ( - "Currently can't handle equality constraints yet.") - repn = self.generate_repn_param(instance, c) - assert repn.is_linear(), ( - "Constraint {} should be linear in " - "unc. parameters".format(c.name)) - - # Generate robust counterpart - det = quicksum(x[0]*x[1].nominal for x in zip(repn.linear_coefs, - repn.linear_vars)) - det += repn.constant - param_var_dict = {id(param): var - for param, var - in zip(repn.linear_vars, repn.linear_coefs)} - # padding = sqrt( var^T * cov^-1 * var ) - padding = quicksum(param_var_dict[id(param[ind_i])] - * uncset.cov[i][j] - * param_var_dict[id(param[ind_j])] - for i, ind_i in enumerate(param) - for j, ind_j in enumerate(param)) - if c.ctype is Constraint: - # For upper bound: det + padding <= b - if c.has_ub(): - name = c.name + '_counterpart_upper' - if root: - expr = det + sqrt(padding) <= c.upper - counterpart = Constraint(expr=expr) - else: - setattr(instance, - c.name + '_padding', - Var(bounds=(0, float('inf')))) - pvar = getattr(instance, c.name + '_padding') - counterpart = Constraint(expr=det + pvar <= c.upper()) - deterministic = Constraint(expr=padding <= pvar**2) - setattr(instance, c.name + '_det_upper', deterministic) - setattr(instance, name, counterpart) - # For lower bound: det - padding >= b - if c.has_lb(): - name = c.name + '_counterpart_lower' - if root: - expr = det - sqrt(padding) >= c.lower - counterpart = Constraint(expr=expr) - else: - setattr(instance, - c.name + '_padding', - Var(bounds=(0, float('inf')))) - pvar = getattr(instance, c.name + '_padding') - counterpart = Constraint(expr=c.lower() <= det - pvar) - deterministic = Constraint(expr=padding <= pvar**2) - setattr(instance, c.name + '_det_lower', deterministic) - setattr(instance, name, counterpart) - else: - # For minimization: min det + padding - # For maximization: max det - padding - sense = c.sense - name = c.name + '_counterpart' - if root: - expr = det + c.sense*sqrt(padding) - counterpart = Objective(expr=expr, sense=sense) - else: - setattr(instance, - c.name + '_padding', - Var(bounds=(0, float('inf')))) - pvar = getattr(instance, c.name + '_padding') - counterpart = Objective(expr=det + sense*pvar, sense=sense) - deterministic = Constraint(expr=padding <= pvar**2) - setattr(instance, - c.name + '_det', - deterministic) - setattr(instance, name, counterpart) - c.deactivate() - - -@TransformationFactory.register('romodel.polyhedral', - doc="Polyhedral Counterpart") -class PolyhedralTransformation(BaseRobustTransformation): - def _apply_to(self, instance, pao=False): - for c in chain(self.get_uncertain_components(instance), - self.get_uncertain_components(instance, - component=Objective)): - # Collect UncParam and UncSet - param = collect_uncparam(c) - uncset = param.uncset - # Check if uncertainty set is empty - assert not uncset.is_empty(), ("{} does not have any " - "constraints.".format(uncset.name)) - # Check if uncertainty set is polyhedral - if not uncset.is_polyhedral(): - continue - - if c.ctype is Constraint: - assert not c.equality, ( - "Currently can't handle equality constraints yet.") - repn = self.generate_repn_param(instance, c) - assert repn.is_linear(), ( - "Constraint {} should be linear in " - "unc. parameters".format(c.name)) - - # Get coefficients - id_coef_dict = {id(repn.linear_vars[i]): - repn.linear_coefs[i] - for i in range(len(repn.linear_vars))} - c_coefs = [id_coef_dict.get(id(i), 0) for i in param.values()] - cons = repn.constant - - # Add dual constraints d^T * v <= b, P^T * v = x - # Constraint - if c.ctype is Constraint: - # LEQ - if c.has_ub(): - # Create linear dual - if pao: - uncset.obj = Objective(expr=c.body, sense=maximize) - dual = create_linear_dual_from(uncset, - unfixed=param.values()) - o_expr = dual.o.expr - del dual.o - dual.o = Constraint(expr=o_expr <= c.upper) - del uncset.obj - else: - dual = self.create_linear_dual(c_coefs, - c.upper - cons, - uncset.mat, - uncset.rhs) - setattr(instance, c.name + '_counterpart_upper', dual) - # GEQ - if c.has_lb(): - # Create linear dual - if pao: - uncset.obj = Objective(expr=c.body, sense=minimize) - dual = create_linear_dual_from(uncset, - unfixed=param.values()) - o_expr = dual.o.expr - del dual.o - dual.o = Constraint(expr=c.lower <= o_expr) - del uncset.obj - else: - dual = self.create_linear_dual([-1*c for c in c_coefs], - -1*(c.lower - cons), - uncset.mat, - uncset.rhs) - setattr(instance, c.name + '_counterpart_lower', dual) - # Objective - else: - setattr(instance, c.name + '_epigraph', Var()) - epigraph = getattr(instance, c.name + '_epigraph') - sense = c.sense - # Create linear dual - if pao: - uncset.obj = Objective(expr=c.expr, sense=-sense) - dual = create_linear_dual_from(uncset, - unfixed=param.values()) - o_expr = dual.o.expr - del dual.o - dual.o = Constraint(expr=sense*o_expr <= sense*epigraph) - del uncset.obj - else: - dual = self.create_linear_dual([sense*c for c in c_coefs], - sense*(epigraph - cons), - uncset.mat, - uncset.rhs) - - setattr(instance, c.name + '_counterpart', dual) - setattr(instance, - c.name + '_new', - Objective(expr=epigraph, sense=sense)) - # Deactivate original constraint - c.deactivate() - - def create_linear_dual(self, c, b, P, d): - ''' - Robust constraint: - c^T*w <= b for all P*w <= d - ''' - blk = Block() - blk.construct() - n, m = len(P), len(P[0]) - # Add dual variables - blk.var = Var(range(n), within=NonNegativeReals) - # Dual objective - blk.obj = Constraint(expr=quicksum(d[j]*blk.var[j] - for j in range(n)) <= b) - # Dual constraints - blk.cons = ConstraintList() - for i in range(m): - lhs = quicksum(P[j][i]*blk.var[j] for j in range(n)) - if lhs.__class__ not in native_numeric_types or lhs != 0: - blk.cons.add(lhs == c[i]) - - return blk - - -@TransformationFactory.register('romodel.generators', - doc=("Replace uncertain constraints by" - " cutting plane generators")) -class GeneratorTransformation(BaseRobustTransformation): - """ Replace all uncertain constraints by RobustConstraint objects. """ - def _apply_to(self, instance): - self._instance = instance - cons = self.get_uncertain_components(instance) - objs = self.get_uncertain_components(instance, component=Objective) - - tdata = instance._transformation_data['romodel.generators'] - tdata.generators = [] - - for c in cons: - generator = RobustConstraint() - setattr(instance, c.name + '_generator', generator) - - generator.build(c.lower, c.body, c.upper) - tdata.generators.append(generator) - - c.deactivate() - - for o in objs: - generator = RobustConstraint() - setattr(instance, o.name + '_epigraph', Var()) - epigraph = getattr(instance, o.name + '_epigraph') - setattr(instance, o.name + '_generator', generator) - - if o.is_minimizing(): - generator.build(None, o.expr - epigraph, 0) - setattr(instance, - o.name + '_new', - Objective(expr=epigraph, sense=minimize)) - else: - generator.build(0, o.expr - epigraph, None) - setattr(instance, - o.name + '_new', - Objective(expr=epigraph, sense=maximize)) - - tdata.generators.append(generator) - - o.deactivate() - - def get_generator(self, c): - pass - - def get_uncset(self, c): - pass - - -@TransformationFactory.register('romodel.nominal', - doc="Transform robust to nominal model.") -class NominalTransformation(BaseRobustTransformation): - def _apply_to(self, instance): - cons = self.get_uncertain_components(instance) - objs = self.get_uncertain_components(instance, component=Objective) - - smap = {} - - for c in cons: - param = collect_uncparam(c) - for i in param: - smap[id(param[i])] = param[i].nominal - body_nominal = replace_expressions(c.body, smap) - c.set_value((c.lower, body_nominal, c.upper)) - - for o in objs: - param = collect_uncparam(o) - for i in param: - smap[id(param[i])] = param[i].nominal - expr_nominal = replace_expressions(o.expr, smap) - o.expr = expr_nominal - - -@TransformationFactory.register('romodel.unknown', - doc="Check for unknown uncertainty sets.") -class UnknownTransformation(BaseRobustTransformation): - def _apply_to(self, instance): - for c in chain(self.get_uncertain_components(instance), - self.get_uncertain_components(instance, - component=Objective)): - # Collect UncParam and UncSet - param = collect_uncparam(c) - uncset = param.uncset - raise RuntimeError("Cannot reformulate UncSet with unknown " - "geometry: {}".format(uncset.name)) - - -@TransformationFactory.register('romodel.warpedgp', - doc="Reformulate warped Gaussian Process set.") -class WGPTransformation(BaseRobustTransformation): - def _apply_to(self, instance, initialize_wolfe=False): - for c in chain(self.get_uncertain_components(instance), - self.get_uncertain_components(instance, - component=Objective)): - if c.ctype is Constraint and c.equality: - raise RuntimeError( - "'UncParam's cannot appear in equality constraints, " - "unless the constraint also contains adjustable " - "variables.") - - # Collect uncertain parameters and uncertainty set - param = collect_uncparam(c) - uncset = param.uncset - - # Only proceed if uncertainty set is WarpedGPSet - if not uncset.__class__ == WarpedGPSet: - continue - from rogp.util.numpy import _pyomo_to_np, _to_np_obj_array - - repn = self.generate_repn_param(instance, c) - - assert repn.is_linear(), ("Only constraints which are linear in " - "the uncertain parameters are valid for " - "uncertainty set WarpedGPSet.") - - # Constraint may contain subset of elements of UncParam - index_set = [p.index() for p in repn.linear_vars] - x = [[xi] for xi in repn.linear_coefs] - x = _to_np_obj_array(x) - - # Set up Block for Wolfe duality constraints - b = Block() - setattr(instance, c.name + '_counterpart', b) - # Set up extra variables - b.y = Var(param.index_set()) - # Set bounds for extra variables based on UncParam bounds - for i in param: - lb, ub = param[i].lb, param[i].ub - b.y[i].setlb(lb) - b.y[i].setub(ub) - if ((lb is not None and lb is not float('-inf')) - and (ub is not None and ub is not float('inf'))): - b.y[i].value = (ub + lb)/2 - y = _pyomo_to_np(b.y, ind=index_set) - # Setup dual vars based on sense - if c.ctype is Constraint: - if c.has_ub() and not c.has_lb(): - b.u = Var(within=NonPositiveReals) - elif not c.has_ub() and c.has_lb(): - b.u = Var(within=NonNegativeReals) - else: - raise RuntimeError( - "Uncertain constraints with 'WarpedGPSet' " - "currently only support either an upper or a " - "lower bound, not both." - ) - elif c.ctype is Objective: - if c.sense is maximize: - b.u = Var(within=NonPositiveReals) - else: - b.u = Var(within=NonNegativeReals) - u = _pyomo_to_np(b.u) - - # Primal constraint - sub_map = {id(param[i]): b.y[i] for i in param} - if c.ctype is Constraint: - e_new = replace_expressions(c.body, substitution_map=sub_map) - b.primal = Constraint(rule=lambda x: (c.lower, e_new, c.upper)) - else: - e_new = replace_expressions(c.expr, substitution_map=sub_map) - b.primal = Objective(expr=e_new, sense=c.sense) - - # Calculate matrices - gp = uncset.gp - var = uncset.var - if type(var) is dict: - z = [var[i] for i in index_set] - z = _to_np_obj_array(z) - else: - var = var[0] - assert var.index_set() == param.index_set(), ( - "Index set of `UncParam` and `var` in `WarpedGPSet` " - "should be the same. Alternatively use " - "var = {index: [list of vars]}" - ) - z = _pyomo_to_np(var, ind=index_set) - - Sig = gp.predict_cov_latent(z) - mu = gp.predict_mu_latent(z) - dHinv = 1/gp.warp_deriv(y) - dHinv = np.diag(dHinv[:, 0]) - hz = gp.warp(y) - - LHS = np.matmul(Sig, dHinv) - LHS = np.matmul(LHS, x) - RHS = LHS - LHS = LHS + 2*u*(hz - mu) - # Add stationarity condition - b.stationarity = ConstraintList() - for lhs in np.nditer(LHS, ['refs_ok']): - b.stationarity.add(lhs.item() == 0) - RHS = np.matmul(dHinv, RHS) - rhs = np.matmul(x.T, RHS)[0, 0] - lhs = 4*u[0, 0]**2*uncset.F - # Set consistent initial value for u (helps convergence) - if initialize_wolfe: - u0 = np.sqrt(rhs()/4/uncset.F) - if u[0, 0].ub == 0: - u[0, 0].value = -u0 - elif u[0, 0].lb == 0: - u[0, 0].value = u0 - # Dual variable constraint - b.dual = Constraint(expr=lhs == rhs) - - c.deactivate() - - -@TransformationFactory.register('romodel.gp', - doc="Reformulate Gaussian Process set.") -class GPTransformation(BaseRobustTransformation): - def _apply_to(self, instance): - for c in chain(self.get_uncertain_components(instance), - self.get_uncertain_components(instance, - component=Objective)): - if c.ctype is Constraint and c.equality: - raise RuntimeError( - "'UncParam's cannot appear in equality constraints, " - "unless the constraint also contains adjustable " - "variables.") - - # Collect uncertain parameters and uncertainty set - param = collect_uncparam(c) - uncset = param.uncset - - # Only proceed if uncertainty set is GPSet - if not uncset.__class__ == GPSet: - continue - from rogp.util.numpy import _pyomo_to_np, _to_np_obj_array - - repn = self.generate_repn_param(instance, c) - - assert repn.is_linear(), ("Only constraints which are linear in " - "the uncertain parameters are valid for " - "uncertainty set 'GPSet'.") - - # Constraint may contain subset of elements of UncParam - index_set = [p.index() for p in repn.linear_vars] - x = [[xi] for xi in repn.linear_coefs] - x = _to_np_obj_array(x) - - # Calculate matrices - gp = uncset.gp - var = uncset.var - if type(var) is dict: - pass - z = [var[i] for i in index_set] - z = _to_np_obj_array(z) - else: - var = var[0] - assert var.index_set() == param.index_set(), ( - "Index set of `UncParam` and `var` in `GPSet` " - "should be the same. Alternatively use " - "var = {index: [list of vars]}" - ) - z = _pyomo_to_np(var, ind=index_set) - - Sig = gp.predict_cov(z) - mu = gp.predict_mu(z) - - nominal = np.matmul(mu.T, x)[0, 0] + repn.constant - - padding = np.matmul(x.T, Sig) - padding = np.matmul(padding, x) - padding = uncset.F*pyomo_sqrt(padding[0, 0]) - - # Counterpart - if c.ctype is Constraint: - if c.has_ub(): - e_new = nominal + padding <= c.upper - c_new = Constraint(expr=e_new) - setattr(instance, c.name + '_counterpart_upper', c_new) - if c.has_lb(): - e_new = nominal - padding >= c.lower - c_new = Constraint(expr=e_new) - setattr(instance, c.name + '_counterpart_lower', c_new) - else: - e_new = nominal + padding*c.sense - c_new = Objective(expr=e_new, sense=c.sense) - setattr(instance, c.name + '_counterpart', c_new) - - c.deactivate() diff --git a/romodel/reformulate/__init__.py b/romodel/reformulate/__init__.py new file mode 100644 index 0000000..6bfb992 --- /dev/null +++ b/romodel/reformulate/__init__.py @@ -0,0 +1,6 @@ +from .base import (BaseRobustTransformation, GeneratorTransformation, + NominalTransformation, UnknownTransformation) +from .ellipsoidal import EllipsoidalTransformation +from .polyhedral import PolyhedralTransformation +from .gp import GPTransformation +from .warpedgp import WGPTransformation diff --git a/romodel/reformulate/base.py b/romodel/reformulate/base.py new file mode 100644 index 0000000..8b95b86 --- /dev/null +++ b/romodel/reformulate/base.py @@ -0,0 +1,188 @@ +from pyomo.environ import (Constraint, + Var, + Objective, + Block, + maximize, + minimize) +from pyomo.core import Transformation, TransformationFactory +from pyomo.repn import generate_standard_repn +from romodel.visitor import _expression_is_uncertain +from romodel.generator import RobustConstraint +from itertools import chain +from romodel.util import collect_uncparam +from pyomo.core.expr.visitor import replace_expressions + + +class BaseRobustTransformation(Transformation): + def __init__(self): + self._fixed_unc_params = [] + self._fixed_components = {} + + def fix_component(self, instance, component=Var): + fixed = [] + for c in instance.component_data_objects(component, active=True): + if not c.is_fixed(): + c.fix() + fixed.append(c) + self._fixed_components[component] = fixed + + def unfix_component(self, component=Var): + for c in self._fixed_components[component]: + c.unfix() + self._fixed_components[component] = None + + def get_uncertain_components(self, instance, component=Constraint): + """ Return all uncertain components of type `component`. """ + comp_list = instance.component_data_objects(component, active=True) + for c in comp_list: + if hasattr(c, 'body'): + # Constraints + expr = c.body + else: + # Objective + expr = c.expr + + if _expression_is_uncertain(expr): + yield c + + def generate_repn_param(self, cdata): + # Get instance + instance = self._instance + # Fix variables + self.fix_component(instance, component=Var) + if hasattr(cdata, 'body'): + expr = cdata.body + else: + expr = cdata.expr + repn = generate_standard_repn(expr, compute_values=False) + self.unfix_component(component=Var) + return repn + + def _apply_to(self, instance, **kwargs): + self._instance = instance + for c in chain(self.get_uncertain_components(instance), + self.get_uncertain_components(instance, + component=Objective)): + # Collect unc. parameter and unc. set + param = collect_uncparam(c) + uncset = param._uncset + assert uncset is not None, ("No uncertainty set provided for " + "uncertain parameter {}." + .format(param.name)) + + # Check if uncertainty set is empty + assert not uncset.is_empty(), ("{} does not have any " + "constraints.".format(uncset.name)) + # Check if reformulation is applicable to constraint & UncSet + if self._check_applicability(uncset): + # Check constraint + if c.ctype is Constraint: + self._check_constraint(c) + else: + self._check_objective(c) + + counterpart = Block() + setattr(instance, c.name + '_counterpart', counterpart) + self._reformulate(c, param, uncset, counterpart, **kwargs) + + c.deactivate() + + def _reformulate(self, c, param, uncset, counterpart, **kwargs): + raise NotImplementedError + + def _check_applicability(self, c): + raise NotImplementedError + + def _check_constraint(self, c): + return True + + def _check_objective(self, c): + return True + + +@TransformationFactory.register('romodel.generators', + doc=("Replace uncertain constraints by" + " cutting plane generators")) +class GeneratorTransformation(BaseRobustTransformation): + """ Replace all uncertain constraints by RobustConstraint objects. """ + def _apply_to(self, instance): + self._instance = instance + cons = self.get_uncertain_components(instance) + objs = self.get_uncertain_components(instance, component=Objective) + + tdata = instance._transformation_data['romodel.generators'] + tdata.generators = [] + + for c in cons: + generator = RobustConstraint() + setattr(instance, c.name + '_generator', generator) + + generator.build(c.lower, c.body, c.upper) + tdata.generators.append(generator) + + c.deactivate() + + for o in objs: + generator = RobustConstraint() + setattr(instance, o.name + '_epigraph', Var()) + epigraph = getattr(instance, o.name + '_epigraph') + setattr(instance, o.name + '_generator', generator) + + if o.is_minimizing(): + generator.build(None, o.expr - epigraph, 0) + setattr(instance, + o.name + '_new', + Objective(expr=epigraph, sense=minimize)) + else: + generator.build(0, o.expr - epigraph, None) + setattr(instance, + o.name + '_new', + Objective(expr=epigraph, sense=maximize)) + + tdata.generators.append(generator) + + o.deactivate() + + def get_generator(self, c): + pass + + def get_uncset(self, c): + pass + + +@TransformationFactory.register('romodel.nominal', + doc="Transform robust to nominal model.") +class NominalTransformation(BaseRobustTransformation): + def _apply_to(self, instance): + cons = self.get_uncertain_components(instance) + objs = self.get_uncertain_components(instance, component=Objective) + + smap = {} + + for c in cons: + param = collect_uncparam(c) + for i in param: + smap[id(param[i])] = param[i].nominal + body_nominal = replace_expressions(c.body, smap) + c.set_value((c.lower, body_nominal, c.upper)) + + for o in objs: + param = collect_uncparam(o) + for i in param: + smap[id(param[i])] = param[i].nominal + expr_nominal = replace_expressions(o.expr, smap) + o.expr = expr_nominal + + +@TransformationFactory.register('romodel.unknown', + doc="Check for unknown uncertainty sets.") +class UnknownTransformation(BaseRobustTransformation): + def _apply_to(self, instance): + for c in chain(self.get_uncertain_components(instance), + self.get_uncertain_components(instance, + component=Objective)): + # Collect UncParam and UncSet + param = collect_uncparam(c) + uncset = param.uncset + raise RuntimeError("Cannot reformulate UncSet with unknown " + "geometry: {}".format(uncset.name)) diff --git a/romodel/reformulate/ellipsoidal.py b/romodel/reformulate/ellipsoidal.py new file mode 100644 index 0000000..82b9fe5 --- /dev/null +++ b/romodel/reformulate/ellipsoidal.py @@ -0,0 +1,150 @@ +from pyomo.environ import Constraint, Var, quicksum, sqrt, Objective, Block +from romodel.reformulate import BaseRobustTransformation +from pyomo.core import TransformationFactory +from pyomo.repn import generate_standard_repn +from romodel.uncset import UncSet, EllipsoidalSet +import numpy as np + + +@TransformationFactory.register('romodel.ellipsoidal', + doc="Ellipsoidal Counterpart") +class EllipsoidalTransformation(BaseRobustTransformation): + def _check_applicability(self, uncset): + """ + Returns `True` if the reformulation is applicable to `uncset` + + uncset: UncSet + + """ + # Check for library set + if uncset.__class__ == EllipsoidalSet: + return True + # Check generic set + elif uncset.__class__ == UncSet: + first_constraint = True + is_ellipsoidal = False + for c in uncset.component_data_objects(Constraint, active=True): + # make sure set has only one constraint + if first_constraint: + first_constraint = False + else: + return False + # Check if constraint is ellipsoidal + repn = generate_standard_repn(c.body) + if not repn.is_quadratic(): + return False + # TODO: assumes implicitly that there is one UncParam per UncSet + param = repn.quadratic_vars[0][0].parent_component() + # Collect covariance matrix and mean + quadratic_coefs = {(id(x[0]), id(x[1])): c for x, c in + zip(repn.quadratic_vars, repn.quadratic_coefs)} + invcov = [[quadratic_coefs.get((id(param[i]), id(param[j])), 0) + for i in param] for j in param] + invcov = np.array(invcov) + invcov = 1/2*(invcov + invcov.T) + eig, _ = np.linalg.eig(invcov) + cov = np.linalg.inv(invcov) + mean = -1/2*np.matmul(cov, np.array(repn.linear_coefs)) + uncset.mean = {x: mean[i] for i, x in enumerate(param)} + uncset.cov = cov.tolist() + uncset.invcov = invcov + # TODO: need to check repn.constant == mean^T * cov * mean? + + is_ellipsoidal = ((c.has_ub() and np.all(eig > 0)) + or (c.has_lb() and np.all(eig < 0))) + + return is_ellipsoidal + + def _check_constraint(self, c): + """ + Raise an error if the constraint is inappropriate for this + reformulation + + c: Constraint + + """ + assert not c.equality, ( + "Currently can't handle equality constraints yet.") + + def _check_objective(self, o): + """ + Raise an error if the objective is inappropriate for this + reformulation + + c: Objective + + """ + pass + + def _reformulate(self, c, param, uncset, counterpart, root=False): + """ + Reformulate an uncertain constraint or objective + + c: Constraint or Objective + param: UncParam + uncset: UncSet + counterpart: Block + + """ + + # Check constraint/objective + repn = self.generate_repn_param(c) + assert repn.is_linear(), ( + "Constraint {} should be linear in " + "unc. parameters".format(c.name)) + + # Generate robust counterpart + det = quicksum(x[0]*x[1].nominal for x in zip(repn.linear_coefs, + repn.linear_vars)) + det += repn.constant + param_var_dict = {id(param): var + for param, var + in zip(repn.linear_vars, repn.linear_coefs)} + # padding = sqrt( var^T * cov^-1 * var ) + padding = quicksum(param_var_dict[id(param[ind_i])] + * uncset.cov[i][j] + * param_var_dict[id(param[ind_j])] + for i, ind_i in enumerate(param) + for j, ind_j in enumerate(param)) + if c.ctype is Constraint: + # For upper bound: det + padding <= b + if c.has_ub(): + counterpart.upper = Block() + if root: + expr = det + sqrt(padding) <= c.upper + robust = Constraint(expr=expr) + else: + counterpart.upper.padding = Var(bounds=(0, float('inf'))) + pvar = counterpart.upper.padding + robust = Constraint(expr=det + pvar <= c.upper()) + deterministic = Constraint(expr=padding <= pvar**2) + counterpart.upper.det = deterministic + counterpart.upper.rob = robust + # For lower bound: det - padding >= b + if c.has_lb(): + counterpart.lower = Block() + if root: + expr = det - sqrt(padding) >= c.lower + robust = Constraint(expr=expr) + else: + counterpart.lower.padding = Var(bounds=(0, float('inf'))) + pvar = counterpart.lower.padding + robust = Constraint(expr=c.lower() <= det - pvar) + deterministic = Constraint(expr=padding <= pvar**2) + counterpart.lower.det = deterministic + counterpart.lower.rob = robust + else: + # For minimization: min det + padding + # For maximization: max det - padding + sense = c.sense + if root: + expr = det + c.sense*sqrt(padding) + robust = Objective(expr=expr, sense=sense) + else: + counterpart.padding = Var(bounds=(0, float('inf'))) + pvar = counterpart.padding + robust = Objective(expr=det + sense*pvar, sense=sense) + deterministic = Constraint(expr=padding <= pvar**2) + counterpart.det = deterministic + counterpart.rob = robust + diff --git a/romodel/reformulate/gp.py b/romodel/reformulate/gp.py new file mode 100644 index 0000000..17b88c4 --- /dev/null +++ b/romodel/reformulate/gp.py @@ -0,0 +1,96 @@ +from pyomo.environ import Constraint, Objective, TransformationFactory +from pyomo.environ import sqrt as pyomo_sqrt +from romodel.uncset import GPSet +from romodel.reformulate import BaseRobustTransformation +import numpy as np + + +@TransformationFactory.register('romodel.gp', + doc="Reformulate Gaussian Process set.") +class GPTransformation(BaseRobustTransformation): + def _check_applicability(self, uncset): + """ + Returns `True` if the reformulation is applicable to `uncset` + + uncset: UncSet + + """ + # Only proceed if uncertainty set is GPSet + return uncset.__class__ == GPSet + + def _check_constraint(self, c): + """ + Raise an error if the constraint is inappropriate for this + reformulation + + c: Constraint + + """ + if c.equality: + raise RuntimeError( + "'UncParam's cannot appear in equality constraints, " + "unless the constraint also contains adjustable " + "variables.") + + def _reformulate(self, c, param, uncset, counterpart): + """ + Reformulate an uncertain constraint or objective + + c: Constraint or Objective + param: UncParam + uncset: UncSet + counterpart: Block + + """ + from rogp.util.numpy import _pyomo_to_np, _to_np_obj_array + + repn = self.generate_repn_param(c) + + assert repn.is_linear(), ("Only constraints which are linear in " + "the uncertain parameters are valid for " + "uncertainty set 'GPSet'.") + + # Constraint may contain subset of elements of UncParam + index_set = [p.index() for p in repn.linear_vars] + x = [[xi] for xi in repn.linear_coefs] + x = _to_np_obj_array(x) + + # Calculate matrices + gp = uncset.gp + var = uncset.var + if type(var) is dict: + pass + z = [var[i] for i in index_set] + z = _to_np_obj_array(z) + else: + var = var[0] + assert var.index_set() == param.index_set(), ( + "Index set of `UncParam` and `var` in `GPSet` " + "should be the same. Alternatively use " + "var = {index: [list of vars]}" + ) + z = _pyomo_to_np(var, ind=index_set) + + Sig = gp.predict_cov(z) + mu = gp.predict_mu(z) + + nominal = np.matmul(mu.T, x)[0, 0] + repn.constant + + padding = np.matmul(x.T, Sig) + padding = np.matmul(padding, x) + padding = uncset.F*pyomo_sqrt(padding[0, 0]) + + # Counterpart + if c.ctype is Constraint: + if c.has_ub(): + e_new = nominal + padding <= c.upper + c_new = Constraint(expr=e_new) + counterpart.upper = c_new + if c.has_lb(): + e_new = nominal - padding >= c.lower + c_new = Constraint(expr=e_new) + counterpart.lower = c_new + else: + e_new = nominal + padding*c.sense + c_new = Objective(expr=e_new, sense=c.sense) + counterpart.obj = c_new diff --git a/romodel/reformulate/polyhedral.py b/romodel/reformulate/polyhedral.py new file mode 100644 index 0000000..be09e15 --- /dev/null +++ b/romodel/reformulate/polyhedral.py @@ -0,0 +1,172 @@ +from pyomo.environ import (Constraint, + Var, + quicksum, + Objective, + maximize, + minimize, + ConstraintList, + NonNegativeReals, + Block, + native_numeric_types) +from pyomo.core import TransformationFactory +from romodel.duality import create_linear_dual_from +from romodel.reformulate import BaseRobustTransformation +from pyomo.repn import generate_standard_repn +from romodel.uncset import UncSet, PolyhedralSet +from romodel.util import collect_uncparam + + +@TransformationFactory.register('romodel.polyhedral', + doc="Polyhedral Counterpart") +class PolyhedralTransformation(BaseRobustTransformation): + def _check_applicability(self, uncset): + """ + Returns `True` if the reformulation is applicable to `uncset` + + uncset: UncSet + + """ + # Check for library set + if uncset.__class__ == PolyhedralSet: + return True + # Check generic set: + elif uncset.__class__ == UncSet: + mat = [] + rhs = [] + for c in uncset.component_data_objects(Constraint, active=True): + # Generate standard repn + repn = generate_standard_repn(c.body) + param = collect_uncparam(c) + # If uncertainty set contains a non-linear constraint it's not + # polyhedral. + if not repn.is_linear(): + return False + coef_dict = {id(x): y for x, y in zip(repn.linear_vars, + repn.linear_coefs)} + if c.has_ub(): + mat.append([coef_dict.get(id(param[i]), 0) for i in param]) + rhs.append(c.upper - repn.constant) + if c.has_lb(): + mat.append([-coef_dict.get(id(param[i]), 0) for i in param]) + rhs.append(repn.constant - c.lower) + + uncset.mat = mat + uncset.rhs = rhs + + return True + + return False + + def _check_constraint(self, c): + """ + Raise an error if the constraint is inappropriate for this + reformulation + + c: Constraint + + """ + assert not c.equality, ( + "Currently can't handle equality constraints yet.") + + def _reformulate(self, c, param, uncset, counterpart, pao=False): + """ + Reformulate an uncertain constraint or objective + + c: Constraint or Objective + param: UncParam + uncset: UncSet + counterpart: Block + + """ + + repn = self.generate_repn_param(c) + assert repn.is_linear(), ( + "Constraint {} should be linear in " + "unc. parameters".format(c.name)) + # Get coefficients + id_coef_dict = {id(repn.linear_vars[i]): + repn.linear_coefs[i] + for i in range(len(repn.linear_vars))} + c_coefs = [id_coef_dict.get(id(i), 0) for i in param.values()] + cons = repn.constant + + # Add dual constraints d^T * v <= b, P^T * v = x + # Constraint + if c.ctype is Constraint: + # LEQ + if c.has_ub(): + # Create linear dual + if pao: + uncset.obj = Objective(expr=c.body, sense=maximize) + dual = create_linear_dual_from(uncset, + unfixed=param.values()) + o_expr = dual.o.expr + del dual.o + dual.o = Constraint(expr=o_expr <= c.upper) + del uncset.obj + else: + dual = self.create_linear_dual(c_coefs, + c.upper - cons, + uncset.mat, + uncset.rhs) + counterpart.upper = dual + # GEQ + if c.has_lb(): + # Create linear dual + if pao: + uncset.obj = Objective(expr=c.body, sense=minimize) + dual = create_linear_dual_from(uncset, + unfixed=param.values()) + o_expr = dual.o.expr + del dual.o + dual.o = Constraint(expr=c.lower <= o_expr) + del uncset.obj + else: + dual = self.create_linear_dual([-1*c for c in c_coefs], + -1*(c.lower - cons), + uncset.mat, + uncset.rhs) + counterpart.lower = dual + # Objective + else: + counterpart.epigraph = Var() + epigraph = counterpart.epigraph + sense = c.sense + # Create linear dual + if pao: + uncset.obj = Objective(expr=c.expr, sense=-sense) + dual = create_linear_dual_from(uncset, + unfixed=param.values()) + o_expr = dual.o.expr + del dual.o + dual.o = Constraint(expr=sense*o_expr <= sense*epigraph) + del uncset.obj + else: + dual = self.create_linear_dual([sense*c for c in c_coefs], + sense*(epigraph - cons), + uncset.mat, + uncset.rhs) + counterpart.dual = dual + counterpart.obj = Objective(expr=epigraph, sense=sense) + + def create_linear_dual(self, c, b, P, d): + ''' + Robust constraint: + c^T*w <= b for all P*w <= d + ''' + blk = Block() + blk.construct() + n, m = len(P), len(P[0]) + # Add dual variables + blk.var = Var(range(n), within=NonNegativeReals) + # Dual objective + blk.obj = Constraint(expr=quicksum(d[j]*blk.var[j] + for j in range(n)) <= b) + # Dual constraints + blk.cons = ConstraintList() + for i in range(m): + lhs = quicksum(P[j][i]*blk.var[j] for j in range(n)) + if lhs.__class__ not in native_numeric_types or lhs != 0: + blk.cons.add(lhs == c[i]) + + return blk diff --git a/romodel/reformulate/warpedgp.py b/romodel/reformulate/warpedgp.py new file mode 100644 index 0000000..f313445 --- /dev/null +++ b/romodel/reformulate/warpedgp.py @@ -0,0 +1,147 @@ +from pyomo.environ import (Constraint, + Var, + Objective, + maximize, + ConstraintList, + NonNegativeReals, + NonPositiveReals) +from pyomo.core import TransformationFactory +from romodel.uncset import WarpedGPSet +from pyomo.core.expr.visitor import replace_expressions +from romodel.reformulate import BaseRobustTransformation +import numpy as np + + +@TransformationFactory.register('romodel.warpedgp', + doc="Reformulate warped Gaussian Process set.") +class WGPTransformation(BaseRobustTransformation): + def _check_applicability(self, uncset): + """ + Returns `True` if the reformulation is applicable to `uncset` + + uncset: UncSet + + """ + return uncset.__class__ == WarpedGPSet + + def _check_constraint(self, c): + """ + Raise an error if the constraint is inappropriate for this + reformulation + + c: Constraint + + """ + if c.equality: + raise RuntimeError( + "'UncParam's cannot appear in equality constraints, " + "unless the constraint also contains adjustable " + "variables.") + + def _reformulate(self, c, param, uncset, counterpart, initialize_wolfe=False): + """ + Reformulate an uncertain constraint or objective + + c: Constraint or Objective + param: UncParam + uncset: UncSet + counterpart: Block + + """ + + # Only proceed if uncertainty set is WarpedGPSet + from rogp.util.numpy import _pyomo_to_np, _to_np_obj_array + + repn = self.generate_repn_param(c) + + assert repn.is_linear(), ("Only constraints which are linear in " + "the uncertain parameters are valid for " + "uncertainty set WarpedGPSet.") + + # Constraint may contain subset of elements of UncParam + index_set = [p.index() for p in repn.linear_vars] + x = [[xi] for xi in repn.linear_coefs] + x = _to_np_obj_array(x) + + # Set up Block for Wolfe duality constraints + b = counterpart + # Set up extra variables + b.y = Var(param.index_set()) + # Set bounds for extra variables based on UncParam bounds + for i in param: + lb, ub = param[i].lb, param[i].ub + b.y[i].setlb(lb) + b.y[i].setub(ub) + if ((lb is not None and lb is not float('-inf')) + and (ub is not None and ub is not float('inf'))): + b.y[i].value = (ub + lb)/2 + y = _pyomo_to_np(b.y, ind=index_set) + # Setup dual vars based on sense + if c.ctype is Constraint: + if c.has_ub() and not c.has_lb(): + b.u = Var(within=NonPositiveReals) + elif not c.has_ub() and c.has_lb(): + b.u = Var(within=NonNegativeReals) + else: + raise RuntimeError( + "Uncertain constraints with 'WarpedGPSet' " + "currently only support either an upper or a " + "lower bound, not both." + ) + elif c.ctype is Objective: + if c.sense is maximize: + b.u = Var(within=NonPositiveReals) + else: + b.u = Var(within=NonNegativeReals) + u = _pyomo_to_np(b.u) + + # Primal constraint + sub_map = {id(param[i]): b.y[i] for i in param} + if c.ctype is Constraint: + e_new = replace_expressions(c.body, substitution_map=sub_map) + b.primal = Constraint(rule=lambda x: (c.lower, e_new, c.upper)) + else: + e_new = replace_expressions(c.expr, substitution_map=sub_map) + b.primal = Objective(expr=e_new, sense=c.sense) + + # Calculate matrices + gp = uncset.gp + var = uncset.var + if type(var) is dict: + z = [var[i] for i in index_set] + z = _to_np_obj_array(z) + else: + var = var[0] + assert var.index_set() == param.index_set(), ( + "Index set of `UncParam` and `var` in `WarpedGPSet` " + "should be the same. Alternatively use " + "var = {index: [list of vars]}" + ) + z = _pyomo_to_np(var, ind=index_set) + + Sig = gp.predict_cov_latent(z) + mu = gp.predict_mu_latent(z) + dHinv = 1/gp.warp_deriv(y) + dHinv = np.diag(dHinv[:, 0]) + hz = gp.warp(y) + + LHS = np.matmul(Sig, dHinv) + LHS = np.matmul(LHS, x) + RHS = LHS + LHS = LHS + 2*u*(hz - mu) + # Add stationarity condition + b.stationarity = ConstraintList() + for lhs in np.nditer(LHS, ['refs_ok']): + b.stationarity.add(lhs.item() == 0) + RHS = np.matmul(dHinv, RHS) + rhs = np.matmul(x.T, RHS)[0, 0] + lhs = 4*u[0, 0]**2*uncset.F + # Set consistent initial value for u (helps convergence) + if initialize_wolfe: + u0 = np.sqrt(rhs()/4/uncset.F) + if u[0, 0].ub == 0: + u[0, 0].value = -u0 + elif u[0, 0].lb == 0: + u[0, 0].value = u0 + # Dual variable constraint + b.dual = Constraint(expr=lhs == rhs) diff --git a/romodel/tests/test_reformulate.py b/romodel/tests/test_reformulate.py index ba4e4a6..776439e 100644 --- a/romodel/tests/test_reformulate.py +++ b/romodel/tests/test_reformulate.py @@ -16,11 +16,11 @@ class TestReformulation(unittest.TestCase): def test_polyhedral(self): m = romodel.examples.Knapsack() m.w.uncset = m.P - self.assertTrue(m.P.is_polyhedral()) - self.assertTrue(m.Plib.is_polyhedral()) - self.assertFalse(m.E.is_polyhedral()) - self.assertFalse(m.Elib.is_polyhedral()) t = PolyhedralTransformation() + self.assertTrue(t._check_applicability(m.P)) + self.assertTrue(t._check_applicability(m.Plib)) + self.assertFalse(t._check_applicability(m.E)) + self.assertFalse(t._check_applicability(m.Elib)) t.apply_to(m) solver = pe.SolverFactory('gurobi_direct') solver.solve(m) @@ -52,7 +52,8 @@ def test_polyhedral_cons_lb(self): t = ro.PolyhedralTransformation() t.apply_to(m) self.assertFalse(m.cons.active) - self.assertTrue(hasattr(m, 'cons_counterpart_lower')) + self.assertTrue(hasattr(m, 'cons_counterpart')) + self.assertTrue(hasattr(m.cons_counterpart, 'lower')) def test_polyhedral_cons_ub(self): m = pe.ConcreteModel() @@ -69,7 +70,8 @@ def test_polyhedral_cons_ub(self): t = ro.PolyhedralTransformation() t.apply_to(m) self.assertFalse(m.cons.active) - self.assertTrue(hasattr(m, 'cons_counterpart_upper')) + self.assertTrue(hasattr(m, 'cons_counterpart')) + self.assertTrue(hasattr(m.cons_counterpart, 'upper')) def test_polyhedral_obj_min(self): m = pe.ConcreteModel() @@ -87,8 +89,9 @@ def test_polyhedral_obj_min(self): t.apply_to(m) self.assertFalse(m.obj.active) self.assertTrue(hasattr(m, 'obj_counterpart')) - self.assertTrue(hasattr(m, 'obj_new')) - self.assertIs(m.obj_new.sense, pe.minimize) + self.assertTrue(hasattr(m.obj_counterpart, 'obj')) + self.assertTrue(hasattr(m.obj_counterpart, 'dual')) + self.assertIs(m.obj_counterpart.obj.sense, pe.minimize) def test_polyhedral_obj_max(self): m = pe.ConcreteModel() @@ -106,19 +109,20 @@ def test_polyhedral_obj_max(self): t.apply_to(m) self.assertFalse(m.obj.active) self.assertTrue(hasattr(m, 'obj_counterpart')) - self.assertTrue(hasattr(m, 'obj_new')) - self.assertIs(m.obj_new.sense, pe.maximize) + self.assertTrue(hasattr(m.obj_counterpart, 'obj')) + self.assertTrue(hasattr(m.obj_counterpart, 'dual')) + self.assertIs(m.obj_counterpart.obj.sense, pe.maximize) @unittest.skipIf('gurobi_direct' not in solvers, 'gurobi_direct not available') def test_ellipsoidal(self): m = romodel.examples.Knapsack() m.w.uncset = m.E - self.assertFalse(m.P.is_ellipsoidal()) - self.assertFalse(m.Plib.is_ellipsoidal()) - self.assertTrue(m.E.is_ellipsoidal()) - self.assertTrue(m.Elib.is_ellipsoidal()) t = EllipsoidalTransformation() + self.assertFalse(t._check_applicability(m.P)) + self.assertFalse(t._check_applicability(m.Plib)) + self.assertTrue(t._check_applicability(m.E)) + self.assertTrue(t._check_applicability(m.Elib)) t.apply_to(m) solver = pe.SolverFactory('gurobi_direct') solver.options['NonConvex'] = 2 @@ -154,7 +158,9 @@ def test_ellipsoidal_cons_lb(self): t = ro.EllipsoidalTransformation() t.apply_to(m, root=False) self.assertFalse(m.cons.active) - self.assertTrue(hasattr(m, 'cons_counterpart_lower')) + self.assertTrue(hasattr(m, 'cons_counterpart')) + self.assertTrue(hasattr(m.cons_counterpart, 'lower')) + self.assertTrue(hasattr(m.cons_counterpart.lower, 'rob')) def test_ellipsoidal_objective_min(self): m = pe.ConcreteModel() @@ -174,8 +180,10 @@ def test_ellipsoidal_objective_min(self): t.apply_to(m, root=False) self.assertFalse(m.obj.active) self.assertTrue(hasattr(m, 'obj_counterpart')) - self.assertTrue(hasattr(m, 'obj_padding')) - self.assertIs(m.obj_counterpart.sense, pe.minimize) + self.assertTrue(hasattr(m.obj_counterpart, 'padding')) + self.assertTrue(hasattr(m.obj_counterpart, 'det')) + self.assertTrue(hasattr(m.obj_counterpart, 'rob')) + self.assertIs(m.obj_counterpart.rob.sense, pe.minimize) def test_ellipsoidal_objective_max(self): m = pe.ConcreteModel() @@ -195,9 +203,10 @@ def test_ellipsoidal_objective_max(self): t.apply_to(m, root=False) self.assertFalse(m.obj.active) self.assertTrue(hasattr(m, 'obj_counterpart')) - self.assertTrue(hasattr(m, 'obj_det')) - self.assertTrue(hasattr(m, 'obj_padding')) - self.assertIs(m.obj_counterpart.sense, pe.maximize) + self.assertTrue(hasattr(m.obj_counterpart, 'padding')) + self.assertTrue(hasattr(m.obj_counterpart, 'det')) + self.assertTrue(hasattr(m.obj_counterpart, 'rob')) + self.assertIs(m.obj_counterpart.rob.sense, pe.maximize) def test_ellipsoidal_cons_lb_root(self): m = pe.ConcreteModel() @@ -216,7 +225,9 @@ def test_ellipsoidal_cons_lb_root(self): t = ro.EllipsoidalTransformation() t.apply_to(m, root=True) self.assertFalse(m.cons.active) - self.assertTrue(hasattr(m, 'cons_counterpart_lower')) + self.assertTrue(hasattr(m, 'cons_counterpart')) + self.assertTrue(hasattr(m.cons_counterpart, 'lower')) + self.assertTrue(hasattr(m.cons_counterpart.lower, 'rob')) def test_ellipsoidal_obj_min_root(self): m = pe.ConcreteModel() @@ -236,9 +247,9 @@ def test_ellipsoidal_obj_min_root(self): t.apply_to(m, root=True) self.assertFalse(m.obj.active) self.assertTrue(hasattr(m, 'obj_counterpart')) - self.assertFalse(hasattr(m, 'obj_padding')) - self.assertFalse(hasattr(m, 'obj_det')) - self.assertIs(m.obj_counterpart.sense, pe.minimize) + self.assertFalse(hasattr(m.obj_counterpart, 'det')) + self.assertTrue(hasattr(m.obj_counterpart, 'rob')) + self.assertIs(m.obj_counterpart.rob.sense, pe.minimize) def test_ellipsoidal_obj_max_root(self): m = pe.ConcreteModel() @@ -258,9 +269,9 @@ def test_ellipsoidal_obj_max_root(self): t.apply_to(m, root=True) self.assertFalse(m.obj.active) self.assertTrue(hasattr(m, 'obj_counterpart')) - self.assertFalse(hasattr(m, 'obj_det')) - self.assertFalse(hasattr(m, 'obj_padding')) - self.assertIs(m.obj_counterpart.sense, pe.maximize) + self.assertFalse(hasattr(m.obj_counterpart, 'det')) + self.assertTrue(hasattr(m.obj_counterpart, 'rob')) + self.assertIs(m.obj_counterpart.rob.sense, pe.maximize) # def test_ellipsoidal_lib_root(self): # m = romodel.examples.Knapsack() diff --git a/romodel/uncset/base.py b/romodel/uncset/base.py index c27626b..33c00e4 100644 --- a/romodel/uncset/base.py +++ b/romodel/uncset/base.py @@ -40,63 +40,6 @@ def __init__(self, *args, **kwargs): self._lib = False - def is_ellipsoidal(self): - # TODO: assumes there is only one constraint on UncSet - for c in self.component_data_objects(Constraint, active=True): - repn = generate_standard_repn(c.body) - if not repn.is_quadratic(): - return False - # TODO: assumes implicitly that there is one UncParam per UncSet - param = repn.quadratic_vars[0][0].parent_component() - # Collect covariance matrix and mean - quadratic_coefs = {(id(x[0]), id(x[1])): c for x, c in - zip(repn.quadratic_vars, repn.quadratic_coefs)} - invcov = [[quadratic_coefs.get((id(param[i]), id(param[j])), 0) - for i in param] for j in param] - invcov = np.array(invcov) - invcov = 1/2*(invcov + invcov.T) - eig, _ = np.linalg.eig(invcov) - cov = np.linalg.inv(invcov) - mean = -1/2*np.matmul(cov, np.array(repn.linear_coefs)) - self.mean = {x: mean[i] for i, x in enumerate(param)} - self.cov = cov.tolist() - self.invcov = invcov - # TODO: need to check repn.constant == mean^T * cov * mean? - - return ((c.has_ub() and np.all(eig > 0)) - or (c.has_lb() and np.all(eig < 0))) - - def is_polyhedral(self): - mat = [] - rhs = [] - param = None - for c in self.component_data_objects(Constraint, active=True): - # Collect uncertain parameter - for p in identify_parent_components(c.body, [UncParam]): - if param is None: - param = p - else: - assert param is p, ("Uncertainty set {} should " - "only contain one UncParam " - "component.".format(self.name)) - # Generate standard repn - repn = generate_standard_repn(c.body) - # If uncertainty set contains a non-linear constraint it's not - # polyhedral. - if not repn.is_linear(): - return False - coef_dict = {id(x): y for x, y in zip(repn.linear_vars, - repn.linear_coefs)} - if c.has_ub(): - mat.append([coef_dict.get(id(param[i]), 0) for i in param]) - rhs.append(c.upper - repn.constant) - if c.has_lb(): - mat.append([-coef_dict.get(id(param[i]), 0) for i in param]) - rhs.append(repn.constant - c.lower) - self.mat = mat - self.rhs = rhs - return True - def is_empty(self): if self._lib: return False diff --git a/romodel/uncset/ellipsoidal.py b/romodel/uncset/ellipsoidal.py index 6a966b5..59377d4 100644 --- a/romodel/uncset/ellipsoidal.py +++ b/romodel/uncset/ellipsoidal.py @@ -25,9 +25,3 @@ def generate_cons_from_lib(self, param): * invcov[i, j] * (param[ind_j] - self.mean[j])) yield None, expr, self.rhs - - def is_ellipsoidal(self): - return True - - def is_polyhedral(self): - return False diff --git a/romodel/uncset/polyhedral.py b/romodel/uncset/polyhedral.py index ed25a47..5a2b4a5 100644 --- a/romodel/uncset/polyhedral.py +++ b/romodel/uncset/polyhedral.py @@ -19,9 +19,3 @@ def generate_cons_from_lib(self, param): quicksum(row[j]*param[ind] for j, ind in enumerate(param)), self.rhs[i]) - - def is_polyhedral(self): - return True - - def is_ellipsoidal(self): - return False