Sympy中的Steinhart-Hart热敏电阻校准

Posted 2021-04-29

基于热敏电阻模型的Steinhart-Hart方程，我为NTC热敏电阻校准提供了简单的三点代码：

from sympy import *

var('R1 R2 R3')
var('T1 T2 T3')
var('A  B  C')

#Ametherm DG103395
T1,R1 = 273.150,31991.6
T2,R2 = 323.150, 3641.0
T3,R3 = 373.150,  686.2

equations = [
    Eq( A + B*ln(R1) + C*ln(R1)**3 , 1.0/T1 ),
    Eq( A + B*ln(R2) + C*ln(R2)**3 , 1.0/T2 ),
    Eq( A + B*ln(R3) + C*ln(R3)**3 , 1.0/T3 ),
]

print solve(equations, (A,B,C))

理想情况下，应该以minimalistic矩阵形式将其推广到多点，例如

T = Matrix([ 273.150, 323.150, 373.150, ....])
R = Matrix([1991.6,  3641.0,   686.2,   ....])
K = MatrixSymbol('K', 3, 1)

print solve( Eq( [1,ln(?),ln(?)**3]*K - 1/T]) )

通过在循环中创建方程式来扩展原始代码非常简单。但它看起来很灰。

是否存在干净且正确的矩阵方式？

UPDATE：看来我自己找到了解决方案。也许可以减少更多？

from sympy import *

#Ametherm DG103395
T = Matrix([  273.150,   298.150, 323.150,  373.150, ])
R = Matrix([31991.6,   10000.0,  3641.0,   686.2,    ])

M = Matrix.vstack( *R.applyfunc( lambda x: Matrix([[1, ln(x), ln(x)**3]]) ) )

ABC = M.solve_least_squares( T.applyfunc(lambda x: 1/x) )

print ABC

答案

这是一种编写方法，不使用矩阵，而是仅使用列表：

from sympy import *

T = [273.150, 323.150, 373.150]
R = [31991.6, 3641.0, 686.2]
n = len(T)
assert n == len(R), "R and T need to have the same length"

A, B, C = symbols("A B C", real=True)

equations = [ Eq( A + B*ln(R[i]) + C*ln(R[i])**3 , 1.0/T[i] ) for i in range(n) ]
print (solve(equations, (A,B,C)))

输出：{A: 0.00115679797363983, B: 0.000227813584600384, C: 1.26349943638314e-7}

请注意，您不需要将T和R声明为SymPy符号。它们是常量，当SymPy在混合公式中遇到它们时，它们会自动转换。