Source code for pyomo.dae.integral

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2022
#  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.deprecation import RenamedClass
from pyomo.core.base.component import ModelComponentFactory
from pyomo.core.base.indexed_component import rule_wrapper
from pyomo.core.base.expression import (Expression,
                                        _GeneralExpressionData,
                                        ScalarExpression,
                                        IndexedExpression)
from pyomo.dae.contset import ContinuousSet
from pyomo.dae.diffvar import DAE_Error

__all__ = ('Integral', )


[docs]@ModelComponentFactory.register("Integral Expression in a DAE model.") class Integral(Expression): """ Represents an integral over a continuous domain The :py:class:`Integral<pyomo.dae.Integral>` component can be used to represent an integral taken over the entire domain of a :py:class:`ContinuousSet<pyomo.dae.ContinuousSet>`. Once every :py:class:`ContinuousSet<pyomo.dae.ContinuousSet>` in a model has been discretized, any integrals in the model will be converted to algebraic equations using the trapezoid rule. Future development will include more sophisticated numerical integration methods. Parameters ---------- *args Every indexing set needed to evaluate the integral expression wrt : :py:class:`ContinuousSet<pyomo.dae.ContinuousSet>` The continuous domain over which the integral is being taken rule : function Function returning the expression being integrated """ def __new__(cls, *args, **kwds): if cls != Integral: return super(Integral, cls).__new__(cls) if len(args) == 0: raise ValueError("Integral must be indexed by a ContinuousSet") elif len(args) == 1: return ScalarIntegral.__new__(ScalarIntegral) else: return IndexedIntegral.__new__(IndexedIntegral) def __init__(self, *args, **kwds): if "wrt" in kwds and "withrespectto" in kwds: raise TypeError( "Cannot specify both 'wrt' and 'withrespectto keywords") wrt = kwds.pop('wrt', None) wrt = kwds.pop('withrespectto', wrt) if wrt is None: # Check to be sure Integral is indexed by single # ContinuousSet and take Integral with respect to that # ContinuousSet if len(args) != 1: raise ValueError( "Integral indexed by multiple ContinuousSets. " "The desired ContinuousSet must be specified using the " "keyword argument 'wrt'") wrt = args[0] if type(wrt) is not ContinuousSet: raise ValueError( "Cannot take the integral with respect to '%s'. Must take an " "integral with respect to a ContinuousSet" % wrt) self._wrt = wrt loc = None for i, s in enumerate(args): if s is wrt: loc = i # Check that the wrt ContinuousSet is in the argument list if loc is None: raise ValueError( "The ContinuousSet '%s' was not found in the indexing sets " "of the Integral" % wrt.name) self.loc = loc # Remove the index that the integral is being expanded over arg = args[0:loc] + args[loc + 1:] # Check that if bounds are given bounds = kwds.pop('bounds', None) if bounds is not None: raise DAE_Error( "Setting bounds on integrals has not yet been implemented. " "Integrals may only be taken over an entire ContinuousSet") # Create integral expression and pass to the expression initialization intexp = kwds.pop('expr', None) intexp = kwds.pop('rule', intexp) if intexp is None: raise ValueError( "Must specify an integral expression") _is_indexed = bool(len(arg)) def _trap_rule(rule, m, *a): ds = sorted(m.find_component(wrt.local_name)) return sum(0.5 * (ds[i + 1] - ds[i]) * (rule(m, * (a[0:loc] + (ds[i + 1],) + a[loc:])) + rule(m, * (a[0:loc] + (ds[i],) + a[loc:]))) for i in range(len(ds) - 1)) # Note that position_map is mapping arguments (block, *args), so # must be 1 more than len(args), and loc has to be offset by one kwds['rule'] = rule_wrapper(intexp, _trap_rule, positional_arg_map=( i for i in range(len(args)+1) if i != loc+1)) kwds.setdefault('ctype', Integral) Expression.__init__(self, *arg, **kwds)
[docs] def get_continuousset(self): """ Return the :py:class:`ContinuousSet<pyomo.dae.ContinuousSet>` the integral is being taken over """ return self._wrt
class ScalarIntegral(ScalarExpression, Integral): """ An integral that will have no indexing sets after applying a numerical integration transformation """ def __init__(self, *args, **kwds): _GeneralExpressionData.__init__(self, None, component=self) Integral.__init__(self, *args, **kwds) def clear(self): self._data = {} def is_fully_discretized(self): """ Checks to see if all ContinuousSets indexing this Integral have been discretized """ if 'scheme' not in self._wrt.get_discretization_info(): return False return True class SimpleIntegral(metaclass=RenamedClass): __renamed__new_class__ = ScalarIntegral __renamed__version__ = '6.0' class IndexedIntegral(IndexedExpression, Integral): """ An integral that will be indexed after applying a numerical integration transformation """ def is_fully_discretized(self): """ Checks to see if all ContinuousSets indexing this Integral have been discretized. """ wrt = self._wrt if 'scheme' not in wrt.get_discretization_info(): return False setlist = [] if self.dim() == 1: setlist = [self.index_set(), ] else: setlist = self.index_set().set_tuple for i in setlist: if i.ctype is ContinuousSet: if 'scheme' not in i.get_discretization_info(): return False return True