Merge pull request #2395 from jsiirola/units-checking-roundoff

Allow relative tolerance when comparing unit dimentionality

Author: jsiirola

pyomo/core/base/units_container.py

@@ -1082,22 +1082,32 @@ def __getattr__(self, item):
     #                                                                  float(conv_offset))
     #     self._pint_registry.define(defn_str)
 
+    def _rel_diff(self, a, b):
+        scale = min(abs(a), abs(b))
+        if scale < 1.:
+            scale = 1.
+        return abs(a - b) / scale
+
     def _equivalent_pint_units(self, a, b, TOL=1e-12):
         if a is b or a == b:
             return True
         base_a = self._pint_registry.get_base_units(a)
         base_b = self._pint_registry.get_base_units(b)
         if base_a[1] != base_b[1]:
-            return False
-        return abs(base_a[0] - base_b[0]) / min(base_a[0], base_b[0]) <= TOL
+            uc_a = base_a[1].dimensionality
+            uc_b = base_b[1].dimensionality
+            for key in uc_a.keys() | uc_b.keys():
+                if self._rel_diff(uc_a.get(key, 0), uc_b.get(key, 0)) >= TOL:
+                    return False
+        return self._rel_diff(base_a[0], base_b[0]) <= TOL
 
     def _equivalent_to_dimensionless(self, a, TOL=1e-12):
         if a is self._pint_dimensionless or a == self._pint_dimensionless:
             return True
         base_a = self._pint_registry.get_base_units(a)
         if not base_a[1].dimensionless:
             return False
-        return abs(base_a[0] - 1.) / min(base_a[0], 1.) <= TOL
+        return self._rel_diff(base_a[0], 1.) <= TOL
 
     def _get_pint_units(self, expr):
         """

pyomo/util/tests/test_check_units.py

@@ -198,5 +198,31 @@ def arcrule(m):
 
         assert_units_consistent(m)
 
+    def test_units_roundoff_error(self):
+        # Issue 2393: this example resulted in roundoff error where the
+        # computed units for var_1 were
+        #(0.9999999999999986,
+        # <Unit('kilogram ** 1 / kelvin / meter ** 1.66533e-16 / second ** 3')>
+        #)
+        m = ConcreteModel()
+        m.var_1 = Var(
+            initialize=400,
+            units=((units.J**0.4) *
+                   (units.kg**0.2) *
+                   (units.W**0.6) /
+                   units.K /
+                   (units.m**2.2) /
+                   (units.Pa**0.2) /
+                   (units.s**0.8))
+        )
+        m.var_1.fix()
+        m.var_2 = Var(
+            initialize=400,
+            units=units.kg/units.s**3/units.K
+        )
+        assert_units_equivalent(m.var_1, m.var_2)
+
+
+
 if __name__ == "__main__":
     unittest.main()