Merge pull request Pyomo#1494 from carldlaird/units-bugfix

Units: Add support for ContinuousSet, clean handling of equality constraints

pyomo/util/check_units.py

@@ -16,11 +16,11 @@
 from pyomo.core.base.units_container import units, UnitsError, UnitExtractionVisitor
 from pyomo.core.base import (Objective, Constraint, Var, Param,
                              Suffix, Set, RangeSet, Block,
-                             ExternalFunction, Expression)
+                             ExternalFunction, Expression,
+                             value)
+from pyomo.dae import ContinuousSet
+from pyomo.mpec import Complementarity
 from pyomo.gdp import Disjunct, Disjunction
-
-from pyomo.gdp import Disjunct
-from pyomo.gdp import Disjunction
 from pyomo.core.expr.template_expr import IndexTemplate
 from pyomo.core.expr.numvalue import native_types
 
@@ -78,18 +78,23 @@ def _assert_units_consistent_constraint_data(condata):
     ConstraintData object are not consistent or are not equivalent
     with each other.
     """
-    if condata.equality:
-        if condata.lower == 0.0:
-            # Pyomo can rearrange expressions, resulting in a value
-            # of 0 for the RHS that does not have units associated
-            # Therefore, if the RHS is 0, we allow it to be unitless
-            # and check the consistency of the body only
-            assert condata.upper == 0.0
-            _assert_units_consistent_expression(condata.body)
-        else:
-            assert_units_equivalent(condata.lower, condata.body)
+    # Pyomo can rearrange expressions, resulting in a value
+    # of 0 for upper or lower that does not have units associated
+    # Therefore, if the lower and/or upper is 0, we allow it to be unitless
+    # and check the consistency of the body only
+    args = list()
+    if condata.lower is not None and value(condata.lower) != 0.0:
+        args.append(condata.lower)
+
+    args.append(condata.body)
+
+    if condata.upper is not None and value(condata.upper) != 0.0:
+        args.append(condata.upper)
+
+    if len(args) == 1:
+        assert_units_consistent(*args)
     else:
-        assert_units_equivalent(condata.lower, condata.body, condata.upper)
+        assert_units_equivalent(*args)
 
 def _assert_units_consistent_property_expr(obj):
     """
@@ -106,6 +111,22 @@ def _assert_units_consistent_expression(expr):
     """
     pyomo_unit, pint_unit = units._get_units_tuple(expr)
 
+# Complementarities that are not in standard form do not
+# current work with the checking code. The Units container
+# should be modified to allow sum and relationals with zero
+# terms (e.g., unitless). Then this code can be enabled.
+#def _assert_units_complementarity(cdata):
+#    """
+#    Raise an exception if any units in either of the complementarity
+#    expressions are inconsistent, and also check the standard block
+#    methods.
+#    """
+#    if cdata._args[0] is not None:
+#        pyomo_unit, pint_unit = units._get_units_tuple(cdata._args[0])
+#    if cdata._args[1] is not None:
+#        pyomo_unit, pint_unit = units._get_units_tuple(cdata._args[1])
+#    _assert_units_consistent_block(cdata)
+
 def _assert_units_consistent_block(obj):
     """
     This method gets all the components from the block
@@ -127,7 +148,11 @@ def _assert_units_consistent_block(obj):
     Disjunct:_assert_units_consistent_block,
     Disjunction: None,
     Block: _assert_units_consistent_block,
-    ExternalFunction: None
+    ExternalFunction: None,
+    ContinuousSet: None, # ToDo: change this when continuous sets have units assigned
+    # complementarities that are not in normal form are not working yet
+    # see comment in test_check_units
+    # Complementarity: _assert_units_complementarity
     }
 
 def assert_units_consistent(obj):
@@ -151,7 +176,11 @@ def assert_units_consistent(obj):
     if objtype in native_types:
         return
     elif obj.is_expression_type() or objtype is IndexTemplate:
-        _assert_units_consistent_expression(obj)
+        try:
+            _assert_units_consistent_expression(obj)
+        except UnitsError:
+            print('Units problem with expression {}'.format(obj))
+            raise
         return
 
     # if object is not in our component handler, raise an exception
@@ -166,6 +195,15 @@ def assert_units_consistent(obj):
     if obj.is_indexed():
         # check all the component data objects
         for cdata in obj.values():
-            handler(cdata)
+            try:
+                handler(cdata)
+            except UnitsError:
+                print('Error in units when checking {}'.format(cdata))
+                raise
     else:
-        handler(obj)
+        try:
+            handler(obj)
+        except UnitsError:
+                print('Error in units when checking {}'.format(obj))
+                raise
+

pyomo/util/tests/test_check_units.py

@@ -13,13 +13,32 @@
 
 import pyutilib.th as unittest
 from pyomo.environ import *
+from pyomo.dae import ContinuousSet
+from pyomo.mpec import Complementarity, complements
+from pyomo.gdp import Disjunct, Disjunction
 from pyomo.core.base.units_container import (
     pint_available, UnitsError,
 )
 from pyomo.util.check_units import assert_units_consistent, assert_units_equivalent, check_units_equivalent
 
+def python_callback_function(arg1, arg2):
+    return 42.0
+
 @unittest.skipIf(not pint_available, 'Testing units requires pint')
 class TestUnitsChecking(unittest.TestCase):
+    def _create_model_and_vars(self):
+        u = units
+        m = ConcreteModel()
+        m.dx = Var(units=u.m, initialize=0.10188943773836046)
+        m.dy = Var(units=u.m, initialize=0.0)
+        m.vx = Var(units=u.m/u.s, initialize=0.7071067769802851)
+        m.vy = Var(units=u.m/u.s, initialize=0.7071067769802851)
+        m.t = Var(units=u.s, bounds=(1e-5,10.0), initialize=0.0024015570927624456)
+        m.theta = Var(bounds=(0, 0.49*3.14), initialize=0.7853981693583533, units=u.radians)
+        m.a = Param(initialize=-32.2, units=u.ft/u.s**2)
+        m.x_unitless = Var()
+        return m
+
     def test_assert_units_consistent_equivalent(self):
         u = units
         m = ConcreteModel()
@@ -134,5 +153,34 @@ def broken(m,i):
         assert_units_consistent(m.vel_con[2]) # check constraint data
         assert_units_consistent(m.unitless_con[2]) # check unitless constraint data
 
+    def test_assert_units_consistent_all_components(self):
+        # test all scalar components consistent
+        u = units
+        m = self._create_model_and_vars()
+        m.obj = Objective(expr=m.dx/m.t - m.vx)
+        m.con = Constraint(expr=m.dx/m.t == m.vx)
+        # vars already added
+        m.exp = Expression(expr=m.dx/m.t - m.vx)
+        m.suff = Suffix(direction=Suffix.LOCAL)
+        # params already added
+        # sets already added
+        m.rs = RangeSet(5)
+        m.disj1 = Disjunct()
+        m.disj1.constraint = Constraint(expr=m.dx/m.t <= m.vx)
+        m.disj2 = Disjunct()
+        m.disj2.constraint = Constraint(expr=m.dx/m.t <= m.vx)
+        m.disjn = Disjunction(expr=[m.disj1, m.disj2])
+        # block tested as part of model
+        m.extfn = ExternalFunction(python_callback_function, units=u.m/u.s, arg_units=[u.m, u.s])
+        m.conext = Constraint(expr=m.extfn(m.dx, m.t) - m.vx==0)
+        m.cset = ContinuousSet(bounds=(0,1))
+
+        # complementarities do not work yet
+        # The expression system removes the u.m since it is multiplied by zero.
+        # We need to change the units_container to allow 0 when comparing units 
+        # m.compl = Complementarity(expr=complements(m.dx/m.t >= m.vx, m.dx == 0*u.m))
+
+        assert_units_consistent(m)
+
 if __name__ == "__main__":
     unittest.main()