如何拟合温度/热曲线的数据？

Posted 2023-03-12

【中文标题】如何拟合温度/热曲线的数据？

【英文标题】：How can fit the data on temperature/thermal profile?

【发布时间】：2019-08-11 17:06:12

【问题描述】：

我有一个包含某个温度曲线的数据集，我想在温度曲线上拟合或映射测量点，如下所示：

停留时间： 30 分钟

斜坡时间： 1 分钟

周期数： 1000 个周期

测点周期：16分钟

测量点可以出现在高区+150或低区-40

注意： T0（初始时间）不明确，因此时间参考不明确，例如。 T0=0。

我已经在 Pandas DataFrame 中获取了数据：

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit

df = pd.read_csv('D:\SOF.csv', header=None)
data = 'A': A[:,0], 'B': B[:,0], 'Temperature': Temperature[:,0],
        'S':S, 'C':C , 'Measurement_Points':MP
dff = pd.DataFrame(data, columns=['A','B','Temperature','S','C','MP'], index = id_set[:,0])
# Temperature's range is [-40,+150]
# MP's range is [0-3000] from 1st MP till last one
MP = int(len(dff)/480) # calculate number of measurement points 
print(MP)
for cycle in range(MP):             
    j = cycle * 480
    #use mean or average of each 480 values from temperature column of DataFrame to pass for fit on Thermal profile
    Mean_temp = np.mean(df['Temperature'].iloc[j:j+480]) # by using Mean from numpy
    #Mean_temp = df.groupby('Temperature').mean() #by using groupby 

到目前为止，我只是根据这个answer和这个post从scipy.optimize找到curve_fit 但我想知道另一方面，拟合过程如何在这里工作我希望温度值 只四舍五入到最接近的 -40 或 +150。 如果有人可以帮助我，我会很好！

更新： 标准周期性热剖面图如下：

预期结果：

更新的数据样本： data

【问题讨论】：

你能提供一些数据吗？！你想要适应的究竟是什么。如果min(T) 和max(T) 是明确的，停留时间和斜坡也是如此，您只需要“拟合”时间偏移，对吧？

你确定坡道上没有点吗？下谷也是30分钟，做一个周期62分钟？

@mikuszefski 正是您在第一条评论中提到的那些功能很清楚，我只需要在热剖面上安装“测量点”，就像您看到的高温红星和低温蓝星一样基于数据，关键是我的数据集中没有直接时间，而是我有多个测量点，每个测量点需要 16 分钟，因此我可以通过 [MP 数量 *((2*30)+( 3 * 1))]。关于您的第二条评论，每个坡道需要 1 分钟，因此根据标准热配置文件，每个周期将是 63 分钟而不是 62 分钟，我更新了图片。

@mikuszefski 另外关于您在注释中提到的第二条评论，因为 T0 不清楚，并且在斜坡持续时间内发生某些测量的可能性很小！我想绘制它们以便我可以跟踪它们或计算它们例如多久发生一次？我想知道某种公式或/和图表最能描述我测量的数据。

似乎并不难。到目前为止，你有什么尝试吗？也将其发布是一个好主意。附言您的第一张图片建议 62。（30 低 +1 上 + 30 高 + 1 下...重复）

【参考方案1】：

这将是我的起点：

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

### to generate test data
def temp( t , low, high, period, ramp ):
    tRed = t % period
    dwell = period / 2. - ramp
    if tRed < dwell:
        out = high
    elif tRed < dwell + ramp:
        out = high - ( tRed - dwell ) / ramp * ( high - low )
    elif tRed < 2 * dwell + ramp:
        out = low
    elif tRed <= period:
        out = low + ( tRed - 2 * dwell - ramp)/ramp * ( high -low )
    else:
        assert 0
    return out + np.random.normal() 

### A continuous function that somewhat fits the data
### but definitively gets the period and levels. 
### The ramp is less well defined
def fit_func( t, low, high, period, s,  delta):
    return  ( high + low ) / 2. + ( high - low )/2. * np.tanh( s * np.sin( 2 * np.pi * ( t - delta ) / period ) )



time1List = np.arange( 300 ) * 16
time2List = np.linspace( 0, 300 * 16, 7213 )
tempList = np.fromiter( ( temp(t - 6.3 , 41, 155, 63.3, 2.05 ) for t in time1List ), np.float )
funcList = np.fromiter( ( fit_func(t , 41, 155, 63.3, 10., 0 ) for t in time2List ), np.float )

sol, err = curve_fit( fit_func, time1List, tempList, [ 40, 150, 63, 10, 0 ] )
print sol

fittedLow, fittedHigh, fittedPeriod, fittedS, fittedOff = sol
realHigh = fit_func( fittedPeriod / 4., *sol)
realLow = fit_func( 3 / 4. * fittedPeriod, *sol)
print "high, low : ", [ realHigh, realLow ]
print "apprx ramp: ", fittedPeriod/( 2 * np.pi * fittedS ) * 2

realAmp = realHigh - realLow
rampX, rampY = zip( *[ [ t, d ] for t, d in zip( time1List, tempList ) if ( ( d < realHigh - 0.05 * realAmp ) and ( d > realLow + 0.05 * realAmp ) ) ] )
topX, topY = zip( *[ [ t, d ] for t, d in zip( time1List, tempList ) if ( ( d > realHigh - 0.05 * realAmp ) ) ] )
botX, botY = zip( *[ [ t, d ] for t, d in zip( time1List, tempList ) if ( ( d < realLow + 0.05 * realAmp ) ) ] )

fig = plt.figure()
ax = fig.add_subplot( 2, 1, 1 )
bx = fig.add_subplot( 2, 1, 2 )

ax.plot( time1List, tempList, marker='x', linestyle='', zorder=100 )
ax.plot( time2List, fit_func( time2List, *sol ), zorder=0 )

bx.plot( time1List, tempList, marker='x', linestyle='' )
bx.plot( time2List, fit_func( time2List, *sol ) )
bx.plot( rampX, rampY, linestyle='', marker='o', markersize=10, fillstyle='none', color='r')
bx.plot( topX, topY, linestyle='', marker='o', markersize=10, fillstyle='none', color='#00FFAA')
bx.plot( botX, botY, linestyle='', marker='o', markersize=10, fillstyle='none', color='#80DD00')
bx.set_xlim( [ 0, 800 ] )
plt.show()

提供：

>> [155.0445024   40.7417905   63.29983807  13.07677546 -26.36945489]
>> high, low :  [155.04450237880076, 40.741790521444436]
>> apprx ramp:  1.540820542195840

有几点需要注意。如果斜坡与停留时间相比较小，我的拟合功能会更好。此外，人们会在这里找到几篇讨论阶跃函数拟合的帖子。通常，由于拟合需要有意义的导数，因此离散函数是一个问题。至少有两种解决方案。 a) 制作连续版本，拟合并根据您的喜好将结果离散化或 b) 提供离散函数和手动连续导数。

编辑

这就是我对您新发布的数据集的处理：

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit, minimize

def partition( inList, n ):
    return zip( *[ iter( inList ) ] * n )

def temp( t, low, high, period, ramp, off ):
    tRed = (t - off ) % period
    dwell = period / 2. - ramp
    if tRed < dwell:
        out = high
    elif tRed < dwell + ramp:
        out = high - ( tRed - dwell ) / ramp * ( high - low )
    elif tRed < 2 * dwell + ramp:
        out = low
    elif tRed <= period:
        out = low + ( tRed - 2 * dwell - ramp)/ramp * ( high -low )
    else:
        assert 0
    return out

def chi2( params, xData=None, yData=None, verbose=False ):
    low, high, period, ramp, off = params
    th = np.fromiter( ( temp( t, low, high, period, ramp, off ) for t in xData ), np.float )
    diff = ( th - yData )
    diff2 = diff**2
    out = np.sum( diff2 )
    if verbose:
        print '-----------'
        print th
        print diff
        print diff2
        print '-----------'
    return out
    # ~ return th

def fit_func( t, low, high, period, s,  delta):
    return  ( high + low ) / 2. + ( high - low )/2. * np.tanh( s * np.sin( 2 * np.pi * ( t - delta ) / period ) )


inData = np.loadtxt('SOF2.csv', skiprows=1, delimiter=',' )
inData2 = inData[ :, 2 ]
xList = np.arange( len(inData2) )
inData480 = partition( inData2, 480 )
xList480 = partition( xList, 480 )
inDataMean = np.fromiter( (np.mean( x ) for x in inData480 ), np.float )
xMean = np.arange( len( inDataMean) ) * 16
time1List = np.linspace( 0, 16 * len(inDataMean), 500 )

sol, err = curve_fit( fit_func, xMean, inDataMean, [ -40, 150, 60, 10, 10 ] )
print sol

# ~ print chi2([-49,155,62.5,1 , 8.6], xMean, inDataMean )
res = minimize( chi2, [-44.12, 150.0, 62.0,  8.015,  12.3 ], args=( xMean, inDataMean ), method='nelder-mead' )
# ~ print res
print res.x

# ~ print chi2( res.x, xMean, inDataMean, verbose=True )
# ~ print chi2( [-44.12, 150.0, 62.0,  8.015,  6.3], xMean, inDataMean, verbose=True )

fig = plt.figure()
ax = fig.add_subplot( 2, 1, 1 )
bx = fig.add_subplot( 2, 1, 2 )

for x,y in zip( xList480, inData480):
    ax.plot( x, y, marker='x', linestyle='', zorder=100 )

bx.plot( xMean, inDataMean , marker='x', linestyle='' )
bx.plot( time1List, fit_func( time1List, *sol ) )

bx.plot( time1List, np.fromiter( ( temp( t , *res.x ) for t in time1List ), np.float) )
bx.plot( time1List, np.fromiter( ( temp( t , -44.12, 150.0, 62.0,  8.015,  12.3 ) for t in time1List ), np.float) )

plt.show()

>> [-49.53569904 166.92138068  62.56131027   1.8547409    8.75673747]
>> [-34.12188737 150.02194584  63.81464913   8.26491754  13.88344623]

如您所见，坡道上的数据点不适合。那么，可能是 16 分钟时间不是那么恒定？这将是一个问题，因为这不是局部 x 错误，而是累积效应。

【讨论】：

首先感谢您的尝试，您说得对，斜坡时间 curve_fit。由于每个测量点有 480 个值，我们可以计算出每 480 个值中的 mean 可以认为是正确的温度，并将它们从 DataFrame 传递给 for-loop ：np.mean(df[' '].iloc[j:j+480]) 或 df.groupby(' ').mean() 基于结构数据样本。我还想知道我们是否可以排除热配置文件功能以进行进一步更改。

@Mario 抱歉，您的示例数据只包含一个跳转，因此没有合适的时间段。取平均值应该不是问题，除非采样时间在斜坡时间范围内，在这种情况下，您可能会隐藏重要信息。不过，我没有得到“热剖面”的东西。什么意思？

@Mario ...但是这不是您需要知道的时期，对吧？你真正想知道的是什么？

我完全同意你的观点，你已经对热曲线进行了编程，但它是基于你生成的数据。我的问题是如何将数据从DataFrame 传递到curve_fit，以便我们可以从原始数据中提取模式。我注意到您已经手动定义了一些值，例如 fit_func time1List tempList [ -40(Low), +150(High), 63(Period), 10(?), 0(?) ] 但是我希望函数从我的数据集中提取这些值。感谢您的详细考虑。

@Mario 热曲线，尤其是我使用的拟合函数是非线性函数，我将用噪声拟合相同的通用数据。原则上curve_fit 只能找到参数，但会卡在局部最小值甚至不收敛。可以通过为参数提供合理的起始值来避免这种情况。这就是我所做的。 10 是描述高原平坦度和斜坡坡度的参数。所以它不是那么重要。 0 只是时间偏移的开始猜测，故意选错了。

【参考方案2】：

如果您只对两个温度级别感兴趣，这可能会有用：

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

inData = np.loadtxt('SOF.csv', skiprows=1, delimiter=',' )

def gauss( x, s ):
    return 1. / np.sqrt( 2. * np.pi * s**2 ) * np.exp( -x**2 / ( 2. * s**2 ) )

def two_peak( x , a1, mu1, s1, a2, mu2, s2 ):
    return a1 * gauss( x - mu1, s1 ) + a2 * gauss( x - mu2, s2 )

fList = inData[ :, 2 ]

nBins = 2 * int( max( fList ) - min( fList ) )
fig = plt.figure()

ax = fig.add_subplot( 2, 1 , 1 )
ax.plot( fList , marker='x' )
bx = fig.add_subplot( 2, 1 , 2 )
histogram, binEdges, _ = bx.hist( fList, bins=nBins )

binCentre = np.fromiter( (  ( a + b ) / 2. for a,b in zip( binEdges[ 1: ], binEdges[ :-1 ] ) ) , np.float )
sol, err = curve_fit( two_peak, binCentre, histogram, [ 120, min( fList ), 1 ] + [ 500, max( fList ), 1 ] )
print sol[1], sol[4]
print sol[2], sol[5]
bx.plot( binCentre, two_peak( binCentre, *sol ) )
bx.set_yscale( 'log' )
bx.set_ylim( [ 1e-0, 5e3] )
plt.show()

提供：

>> -46.01513424923528 150.06381412858244
>> 1.8737971845243133 0.6964990809008554

和

有趣的是，您的非高原数据几乎为零，因此这可能不是由于斜坡，而是不同的影响。