# ___________________________________________________________________________ # # Pyomo: Python Optimization Modeling Objects # Copyright (c) 2008-2024 # National Technology and Engineering Solutions of Sandia, LLC # Under the terms of Contract DE-NA0003525 with National Technology and # Engineering Solutions of Sandia, LLC, the U.S. Government retains certain # rights in this software. # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ from pyomo.common.collections import ComponentSet from pyomo.core.base import Constraint, Block, value from pyomo.dae.set_utils import ( is_explicitly_indexed_by, get_index_set_except, is_in_block_indexed_by, deactivate_model_at, index_warning, ) def get_inconsistent_initial_conditions( model, time, tol=1e-8, t0=None, allow_skip=True, suppress_warnings=False ): """Finds constraints of the model that are implicitly or explicitly indexed by time and checks if they are consistent to within a tolerance at the initial value of time. Args: model: Model whose constraints to check time: Set whose initial condition will be checked tol: Maximum constraint violation t0: Point in time at which to check constraints Returns: List of constraint data objects that were found to be inconsistent. """ if t0 is None: t0 = time.first() inconsistent = ComponentSet() for con in model.component_objects(Constraint, active=True): if not is_explicitly_indexed_by(con, time): continue if is_in_block_indexed_by(con, time): continue info = get_index_set_except(con, time) non_time_set = info['set_except'] index_getter = info['index_getter'] for non_time_index in non_time_set: index = index_getter(non_time_index, t0) try: condata = con[index] except KeyError: # To allow Constraint.Skip if not suppress_warnings: print(index_warning(con.name, index)) if not allow_skip: raise continue if ( value(condata.body) - value(condata.upper) > tol or value(condata.lower) - value(condata.body) > tol ): inconsistent.add(condata) for blk in model.component_objects(Block, active=True): # What if there are time-indexed blocks at multiple levels # of a hierarchy? # My preferred convention is to only check the first (highest- # level) time index, but distinguishing between different-level # time indices is an expensive operation. if not is_explicitly_indexed_by(blk, time): continue if is_in_block_indexed_by(blk, time): continue info = get_index_set_except(blk, time) non_time_set = info['set_except'] index_getter = info['index_getter'] for non_time_index in non_time_set: index = index_getter(non_time_index, t0) blkdata = blk[index] for condata in blkdata.component_data_objects(Constraint, active=True): if ( value(condata.body) - value(condata.upper) > tol or value(condata.lower) - value(condata.body) > tol ): if condata in inconsistent: raise ValueError( '%s has already been visited. The only way this ' 'should happen is if the model has nested time-' 'indexed blocks, which is not supported.' ) inconsistent.add(condata) return list(inconsistent) def solve_consistent_initial_conditions( model, time, solver, tee=False, allow_skip=True, suppress_warnings=False ): """ Solves a model with all Constraints and Blocks deactivated except at the initial value of the Set time. Reactivates Constraints and Blocks that got deactivated. Args: model: Model that will be solved time: Set whose initial conditions will remain active for solve solver: Something that implements a solve method that accepts a model and tee keyword as arguments tee: tee argument that will be sent to solver's solve method allow_skip: If True, KeyErrors due to Constraint.Skip being used will be ignored suppress_warnings: If True, warnings due to ignored KeyErrors will be suppressed Returns: The object returned by the solver's solve method """ # Need to deactivate discretization equations, wrt time, at t == 0 # This is challenging as the only way (to my knowledge) to do this # is to identify_variables in the expression, find the (assume only one?) # DerivativeVar, and access its get_continuousset_list # I would like a get_continuousset_list for discretization equations. # Possibly as a ComponentMap, possibly as an attribute of some new # DiscEquation subclass of Constraint # Until I have this, this function will only work for backward # discretization schemes # Also, would like to be able to check for zero degrees of freedom here scheme = time.get_discretization_info()['scheme'] if scheme != 'LAGRANGE-RADAU' and scheme != 'BACKWARD Difference': raise NotImplementedError('%s discretization scheme is not supported' % scheme) timelist = list(time)[1:] deactivated_dict = deactivate_model_at( model, time, timelist, allow_skip=allow_skip, suppress_warnings=suppress_warnings, ) result = solver.solve(model, tee=tee) for t in timelist: for comp in deactivated_dict[t]: comp.activate() return result