模型使用牛頓冷卻模型

畫出外界溫度曲線

兩端溫度不確定，假定為

其中h等于外界溫度25°，c根據(jù)溫度計算即可

def T(t,v,T15,T6,T7,T89):a=[0,25,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,25]a=[i/v*60 for i in a]m=0a=[m:=m+i for i in a]Tair=25lmd0=-0.5c0=(T15-Tair)/(math.e**(-lmd0*25/v*60))lmd=0.015c=(T89-Tair)/(math.e**(-lmd*339.5/v*60))if t<a[1]:return c0*math.e**(-lmd0*t)+Tairelif t<a[10]:return T15elif t<a[11]:return T15+(T6-T15)/(a[11]-a[10])*(t-a[10])elif t<a[12]:return T6elif t<a[13]:return T6+(T7-T6)/(a[13]-a[12])*(t-a[12])elif t<a[14]:return T7elif t<a[15]:return T7+(T89-T7)/(a[15]-a[14])*(t-a[14])elif t<a[18]:return T89else:return c*math.e**(-lmd*t)+Tairaa=np.linspace(0,400,4001).tolist() plt.plot(aa,[T(i,70,175,195,235,255)for i in aa])

查看題目中給出的數(shù)據(jù)

def shuju():a=[float(i[0])for i in csv.reader(open('./附件.csv'))]b=[float(i[1])for i in csv.reader(open('./附件.csv'))]plt.plot(a,b)plt.show()

可以看到左端數(shù)據(jù)量過小，并且環(huán)境溫度曲線也是假設(shè)的，這個值最后根據(jù)題目數(shù)據(jù)估計，使用pytorch求解會使全局誤差過大或者溫度變化不符合客觀規(guī)律（k為負）

分段求解k值

def K(start,last):a=[float(i[0])for i in csv.reader(open('./附件.csv')) if float(i[0])<last and float(i[0])>start]b=[float(i[1])for i in csv.reader(open('./附件.csv')) if float(i[0])<last and float(i[0])>start]# k=torch.rand([1,1],requires_grad=True)k=torch.tensor([0.1],requires_grad=True)learning_rate=0.01#學習率learning_times=100#循環(huán)次數(shù)criterion=torch.nn.MSELoss()#損失函數(shù)x=torch.tensor(a)#學習for i in range(learning_times):#print(k.data)y_pre=torch.zeros(len(a))#初始化y_prey_true=torch.tensor(b)#這兩行必須放在循環(huán)內(nèi)，不然每次循環(huán)釋放掉了沒法計算梯度for j in range(len(a)):if j==0:y_pre[j]=y_true[j]continuey_pre[j]=y_pre[j-1]+(T(x[j],70,175,195,235,255)-y_pre[j-1])*k*(x[j]-x[j-1])#即牛頓冷卻模型，遞歸求解出所有溫度的預(yù)測值loss=criterion(y_pre,y_true)#計算損失值# loss=(y_true - y_pre).pow(2).mean()if k.grad is not None:k.grad.data.zero_()#重置梯度loss.backward()#反向傳播k.grad/=len(a)**2#上面循環(huán)計算y_pre時每次循環(huán)都積累了一次梯度k.data=k.data-learning_rate*k.gradif i==learning_times-1:last_y_pre=y_pre.tolist()#plt.plot(a,b)#plt.plot(a,last_y_pre)#plt.show()return k.dataprint(K(t_list[1],t_list[11])) print(K(t_list[11],t_list[13])) print(K(t_list[13],t_list[15])) print(K(t_list[15],t_list[18])) print(K(t_list[18],t_list[-1])) # tensor([0.0177]) # tensor([0.0166]) # tensor([0.0243]) # tensor([0.0202]) # tensor([0.0276])

多個局部最優(yōu)不一定是全局最優(yōu)，后面直接將所有數(shù)據(jù)一起計算loss值得出全局最優(yōu)

求解全局最優(yōu)k值

最左端數(shù)據(jù)太少，會受到全局的影響使得k值為負數(shù)，故暫時先不求解最左端的k值，最后再單獨求

def _loss_(start,k1,k2,k3,k4,k5,time,t_list):a=[float(i[0])for i in csv.reader(open('./附件.csv'))if float(i[0])>start]b=[float(i[1])for i in csv.reader(open('./附件.csv'))if float(i[0])>start]x=torch.tensor(a)y_pre=torch.zeros(len(a))#初始化y_prey_true=torch.tensor(b)#這兩行必須放在循環(huán)內(nèi)，不然每次循環(huán)釋放掉了沒法計算梯度for j in range(len(a)):if j==0:y_pre[j]=y_true[j]continueif(x[j]<t_list[1]):continueelif(x[j]<t_list[11]):y_pre[j]=y_pre[j-1]+(T(x[j],70,175,195,235,255)-y_pre[j-1])*k1*(x[j]-x[j-1])elif(x[j]<t_list[13]):y_pre[j]=y_pre[j-1]+(T(x[j],70,175,195,235,255)-y_pre[j-1])*k2*(x[j]-x[j-1])elif(x[j]<t_list[15]):y_pre[j]=y_pre[j-1]+(T(x[j],70,175,195,235,255)-y_pre[j-1])*k3*(x[j]-x[j-1])elif(x[j]<t_list[18]):y_pre[j]=y_pre[j-1]+(T(x[j],70,175,195,235,255)-y_pre[j-1])*k4*(x[j]-x[j-1])else:y_pre[j]=y_pre[j-1]+(T(x[j],70,175,195,235,255)-y_pre[j-1])*k5*(x[j]-x[j-1])#if not time%99:#plt.plot(x.tolist(),y_pre.tolist())#plt.plot(a,b)#plt.show()return torch.nn.MSELoss()(y_pre,y_true)#計算損失值def _k_together():learning_rate=0.01learning_times=1000t_list=[0.0, 21.42857142857143, 47.57142857142857, 51.857142857142854, 78.0, 82.28571428571429, 108.42857142857143, 112.71428571428572, 138.85714285714286, 143.14285714285714, 169.28571428571428,173.57142857142856, 199.7142857142857, 203.99999999999997, 230.1428571428571, 234.4285714285714, 260.57142857142856, 264.85714285714283, 291.0, 295.2857142857143, 321.42857142857144, 325.7142857142857, 351.8571428571429, 373.28571428571433]#t_list是每個溫區(qū)的起始位置和結(jié)束位置除以速度70k1=torch.tensor([0.2],requires_grad=True)k2=torch.tensor([0.2],requires_grad=True)k3=torch.tensor([0.2],requires_grad=True)k4=torch.tensor([0.2],requires_grad=True)k5=torch.tensor([0.2],requires_grad=True)#初始值定位0.2，減少學習次數(shù)for i in range(learning_times):loss=_loss_(t_list[1],k1,k2,k3,k4,k5,i,t_list)if k1.grad is not None:k1.grad.data.zero_()k2.grad.data.zero_()k3.grad.data.zero_()k4.grad.data.zero_()k5.grad.data.zero_()loss.backward()k1.grad/=200**2k2.grad/=70**2k3.grad/=70**2k4.grad/=70**2k5.grad/=100**2#每段中數(shù)據(jù)個數(shù)的平方，這里我估計了一下，并不是準確值k1.data=k1.data-k1.grad*learning_ratek2.data=k2.data-k2.grad*learning_ratek3.data=k3.data-k3.grad*learning_ratek4.data=k4.data-k4.grad*learning_ratek5.data=k5.data-k5.grad*learning_rate#print(k1,k2,k3,k4,k5)_k_together() #tensor([0.0177], requires_grad=True) tensor([0.0171], requires_grad=True) tensor([0.0242], requires_grad=True) tensor([0.0216], requires_grad=True) tensor([0.0276], requires_grad=True) #0.0177 #0.0171 #0.0242 #0.0216 #0.0276

最后單獨計算左端的k值為0.0128（·這個不是很重要，范圍小數(shù)據(jù)少，對整體的影響不大）

分段	爐前	1-5	6	7	8-9	10-爐后
k	0.0128	0.0177	0.0171	0.0242	0.0216	0.0276

第一問的求解

a=[0,25,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,5,30.5,25] a=[i/78*60 for i in a] m=0 a=[m:=m+i for i in a]h=25 lmd0=-0.5 c0=(173-h)/(math.e**(-lmd0*25/78*60)) lmd=0.015 c=(257-h)/(math.e**(-lmd*339.5/78*60)) def T(t):if t<a[1]:return c0*math.e**(-lmd0*t)+helif t<a[10]:return 173elif t<a[11]:return 173+(198-173)/(a[11]-a[10])*(t-a[10])elif t<a[12]:return 198elif t<a[13]:return 198+(230-198)/(a[13]-a[12])*(t-a[12])elif t<a[14]:# return 235return 230elif t<a[15]:return 230+(257-230)/(a[15]-a[14])*(t-a[14])elif t<a[18]:return 257else:return c*math.e**(-lmd*t)+h def k(t):if t<a[1]:return 0.0128elif t<a[11]:return 0.0177elif t<a[13]:return 0.0171elif t<a[15]:return 0.0242elif t<a[18]:return 0.0216else:return 0.0276aa=np.linspace(0,400,4001).tolist() plt.plot(aa,[T(i)for i in aa]) bb=[25] for i in aa:bb.append(bb[-1]+k(i)*(T(i)-bb[-1])*0.1) bb.pop() plt.plot(aa,bb) plt.show()aa=[round(i,1) for i in aa] print(bb[aa.index(round((a[6]+a[5])/2,1))]) print(bb[aa.index(round((a[11]+a[12])/2,1))]) print(bb[aa.index(round((a[13]+a[14])/2,1))]) print(bb[aa.index(round(a[16],1))]) # 128.42713353818752 # 167.691867327485 # 189.18635779121473 # 224.58207564933335

位置	小溫區(qū)3中點	小溫區(qū)6中點	小溫區(qū)7中點	小溫區(qū)8結(jié)束點
溫度	128.42713353818752	167.691867327485	189.18635779121473	224.58207564933335

總結(jié)

以上是生活随笔為你收集整理的2020A炉温曲线第一题pytorch求解（附代码）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。