Verify correctness of Modelica smoothOrder annotation

Question

This is a follow-up question to my previous question: Modelica smoothOrder of cubic polynomial

I am contributing to a large existing library with maybe 50 smoothOrder annotations in functions. These are often piecewise defined function, with different function definitions in the different branches of an if-else statement.

A manual approach for checking the smoothness can be to evaluate and plot the function and all (partial) derivatives close to the switching condition of the if-else, up to the order that is defined. If the values and plots of derivatives are continuous, the annotation is correct.

But this manual approach is quite time consuming, so it would be great to have some kind of automated check. Does this exist, or can anyone share a sample script or model to help getting started?

Answer 1

I have used sympy to do more or less what you described.

The jupyter notebook below computes the partial derivatives of all input equations with respect to the given variable, up to the given order. Then it evaluates the derivatives for a specific value (the switching point in your if conditions). If the output values are equal, the equations should be smooth for that derivative.

# In[1]:

from IPython.display import display
from sympy import init_printing, symbols, diff, simplify
init_printing()


# In[2]:

def der_at(eqs, wrt, nder, at):
    """ Display derivatives of equations up to given order and evaluate for given value """

    ordinal = lambda n: "%d%s" % (n,"tsnrhtdd"[(n//10%10!=1)*(n%10<4)*n%10::4])

    for i in range(0, nder+1):
            
        der_eqs = [diff(eq, wrt, i) for eq in eqs]
        
        der_eqs_subs = [simplify(der_eq.subs([(wrt, at)])) for der_eq in der_eqs]

        print(f"{ordinal(i)}-order derivatives:")
        display(der_eqs)

        print(f"Derivatives evaluated for {wrt}={at}")
        display(der_eqs_subs)


# In[3]:

a, b, c, d, x = symbols('a b c d x')

eq1 = a*x**2 + b*x
eq2 = c*x**3 + d*x**2

der_at(eqs=(eq1, eq2), wrt=x, nder=3, at=0)

After running the code above, you should get the output shown in the image. For the exemplary equations I used, you can see that the we get the same result when we evaluate the original equations for x=0, but not for the first derivatives. Hence the smooth order is 0.