Units not consistent when differentiate returns zero

[dallan-keylogic]

Summary

When trying to use Pyomo's AD to avoid having the user provide derivative expressions, I discovered that, if the resulting derivative is zero, its units get lost and result in dimensional inconsistencies.

It might make sense to treat zero as having units compatible with everything, but I don't know whether that's possible.

Is there a workaround where I can check if the expression returned by differentiate is zero and provide appropriate units manually?

Steps to reproduce the issue

# diff_units.py
import pyomo.environ as pyo
from pyomo.core.expr.calculus.derivatives import Modes, differentiate
from pyomo.util.check_units import assert_units_equivalent

m = pyo.ConcreteModel()

m.temperature = pyo.Var(units=pyo.units.K, initialize=0)
m.density_param = pyo.Param(
    mutable=True,
    initialize=1000,
    units=pyo.units.kg / pyo.units.m ** 3
)

@m.Expression()
def density(blk):
    return blk.density_param


diff_expr = differentiate(
                expr=m.density,
                wrt=m.temperature,
                mode=Modes.reverse_symbolic
            )
print(diff_expr)
assert_units_equivalent(
    diff_expr,
    pyo.units.kg / pyo.units.m ** 3 / pyo.units.K
)

Error Message

Units between dimensionless and kilogram / kelvin / meter ** 3 are not consistent.
  File "C:\Users\[REDACTED]\Desktop\diff_units.py", line 25, in <module>
    assert_units_equivalent(
pyomo.core.base.units_container.UnitsError: Units between dimensionless and kilogram / kelvin / meter ** 3 are not consistent.

Information on your system

Pyomo version: 6.7.2a0 Python version: 3.10.14 Operating system: Windows 10 How Pyomo was installed (PyPI, conda, source): IDAES advanced user installation Solver (if applicable): n/a

[michaelbynum]

Interesting. I've never even considered what happens to units when using differentiate.

@dallan-keylogic - can you provide a little more context about what you are trying to do? Where are these derivatives getting used?

[dallan-keylogic]

@michaelbynum , they're getting used as part of building up an expression for excess molar enthalpy in an activity coefficient model. Initially, I was using AD to build up the entire expression, but it generally builds expressions that are considerably larger than those I can generate by hand. Right now, I'm using it to get derivatives of other user-defined properties to avoid requiring the user to provide the derivative in addition to the property.

v = pyunits.convert(
    getattr(b, pname + "_vol_mol_solvent"), pyunits.m**3 / pyunits.mol
)
dv_dT = differentiate(
    v,
    b.temperature,
    mode=Modes.reverse_symbolic
)

eps = getattr(b, pname + "_relative_permittivity_solvent")
d_eps_dT = getattr(b, pname + "_d_relative_permittivity_solvent_dT")
A_DH = getattr(b, pname + "_A_DH")
Ix = getattr(b, pname + "_ionic_strength")
rho = ClosestApproach
# Temperature derivative of Debye Huckel term
dA_dT = -A_DH / 2 * (
    1/v * dv_dT + 3 * (1/eps) * d_eps_dT
)

[michaelbynum]

Could you use numeric differentiation instead of symbolic? Or you need symbolic?

[dallan-keylogic]

This is in the context of building up a big Pyomo expression, so it has to be symbolic.

[michaelbynum]

I think we could add some logic to the differentiate function that checks if an int or a float is being returned and, if so, makes it a Pyomo expression with units that match what is expected. Thoughts, @jsiirola?

[michaelbynum]

Oh, this is actually because of the way we initialize the derivative map here: https://github.com/Pyomo/pyomo/blob/6365a2526b8d5f067fa0da1e876638c7a7aba38e/pyomo/core/expr/calculus/diff_with_pyomo.py#L423. I think I just need to update that to include units...

[michaelbynum]

No, I think this is probably easier to handle inside of differentiate.

[dallan-keylogic]

Somewhat tangential to this issue, I have another entry in the "zero should be treated as being equivalent with any set of units" file:

pyunits.convert(
      N_A * (q_e)**2/(8 * pi * eps_vacuum * eps_solvent)
      * (1 + b.temperature / eps_solvent * d_eps_solvent_dT)
      * sum(
          b.mole_frac_phase_comp_true[pname, s]
          * abs(cobj(b, s).config.charge) ** 2
          / cobj(b, s).born_radius
          for s in b.params.ion_set
      ),
      to_units=units.ENERGY/units.AMOUNT
)

I have this expression as part of an eNRTL expression, and have just debugged an issue arising when b.params.ion_set is empty. The summation, instead of picking up units of 1/length, returns the integer zero instead, which leads to dimensional inconsistency when the unit conversion is performed.

There's a straightforward solution in this case, however: return zero with correct units if len(b.params.ion_set) == 0, otherwise return this expression.