验证 Modelica smoothOrder 注解的正确性
Verify correctness of Modelica smoothOrder annotation
这是我之前问题的后续问题:
我正在为一个大型现有库做贡献,函数中可能有 50 个 smoothOrder 注释。这些通常是分段定义的函数,在 if-else 语句的不同分支中具有不同的函数定义。
检查平滑度的手动方法可以是评估和绘制函数和所有(偏)导数接近 if-else 的切换条件,直到定义的顺序。如果导数的值和图是连续的,则注释是正确的。
但是这种手动方法非常耗时,所以最好能有某种自动检查。这是否存在,或者任何人都可以共享示例脚本或模型来帮助入门吗?
我使用 sympy 或多或少做了你描述的事情。
下面的 jupyter notebook 计算所有输入方程关于给定变量的偏导数,直到给定阶数。然后它评估特定值的导数(if 条件中的切换点)。如果输出值相等,则方程对于该导数应该是平滑的。
# 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)
在 运行 上面的代码之后,您应该得到如图所示的输出。对于我使用的示例方程,您可以看到当我们计算 x=0
的原始方程时我们得到相同的结果,但不是一阶导数。因此平滑顺序为0.
这是我之前问题的后续问题:
我正在为一个大型现有库做贡献,函数中可能有 50 个 smoothOrder 注释。这些通常是分段定义的函数,在 if-else 语句的不同分支中具有不同的函数定义。
检查平滑度的手动方法可以是评估和绘制函数和所有(偏)导数接近 if-else 的切换条件,直到定义的顺序。如果导数的值和图是连续的,则注释是正确的。
但是这种手动方法非常耗时,所以最好能有某种自动检查。这是否存在,或者任何人都可以共享示例脚本或模型来帮助入门吗?
我使用 sympy 或多或少做了你描述的事情。
下面的 jupyter notebook 计算所有输入方程关于给定变量的偏导数,直到给定阶数。然后它评估特定值的导数(if 条件中的切换点)。如果输出值相等,则方程对于该导数应该是平滑的。
# 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)
在 运行 上面的代码之后,您应该得到如图所示的输出。对于我使用的示例方程,您可以看到当我们计算 x=0
的原始方程时我们得到相同的结果,但不是一阶导数。因此平滑顺序为0.