用GEKKO（二阶微分方程）估计参数

2024-05-23 19:09:44 发布

您现在位置：Python中文网/ 问答频道 /正文
9304 3

我使用的是二阶并行模型：

enter image description here

它有四个参数（初始值Cf0和Cs0、kf和ks），可使用GEKKO和六个实验进行估计。然而，我并没有得到一个很好的拟合，只有实验1和4靠近实验数据。也许我的代码有错误。尽管如此，获得更多控制边界和估算方法信息的任何指导线或链接都将非常有用

enter image description here

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt 
import math as math
import pandas as pd

tm1 = [0, 0.0667,0.5,1,4]
mca1 = [5.68, 3.48, 3.24, 3.36, 2.96]

tm2 = [0, 0.08333,0.5,1,4.25]
mca2 = [5.68, 4.20, 4.04, 4.00, 3.76]

tm3 = [0,0.08333,0.5,1,4.33]
mca3 = [5.68, 4.64, 4.52, 4.56, 4.24]

tm4 = [0,0.08333,0.5,1,4.0833]
mca4 =[18.90,15.4,14.3,15.10,13.50]

tm5 = [0,0.08333,0.5,1,4.5]
mca5 =[18.90, 15.5, 16.30, 16, 14.70]

tm6 = [0,0.08333,0.5,1,4.6667]
mca6 = [18.90, 15.8, 11.70,15.5,12]

df1=pd.DataFrame({'time':tm1,'x1':mca1})
df2=pd.DataFrame({'time':tm2,'x2':mca2})
df3=pd.DataFrame({'time':tm3,'x3':mca3})
df4=pd.DataFrame({'time':tm4,'x4':mca4})
df5=pd.DataFrame({'time':tm5,'x5':mca5})
df6=pd.DataFrame({'time':tm6,'x6':mca6})
df1.set_index('time',inplace=True)
df2.set_index('time',inplace=True)
df3.set_index('time',inplace=True)
df4.set_index('time',inplace=True)
df5.set_index('time',inplace=True)
df6.set_index('time',inplace=True)
#simulation time points
dfx = pd.DataFrame({'time':np.linspace(0,5,101)})
dfx.set_index('time',inplace=True)
#merge dataframes
dfxx=dfx.join(df1,how='outer')
dfxxx=dfxx.join(df2,how='outer')
dfxxxx=dfxxx.join(df3,how='outer')
dfxxxxx=dfxxxx.join(df4,how='outer')
dfxxxxxx=dfxxxxx.join(df5,how='outer')
df=dfxxxxxx.join(df6,how='outer')
# get True (1) or False (0) for measurement
df['meas1']=(df['x1'].values==df['x1'].values).astype(int)
df['meas2']=(df['x2'].values==df['x2'].values).astype(int)
df['meas3']=(df['x3'].values==df['x3'].values).astype(int)
df['meas4']=(df['x4'].values==df['x4'].values).astype(int)
df['meas5']=(df['x5'].values==df['x5'].values).astype(int)
df['meas6']=(df['x6'].values==df['x6'].values).astype(int)
#replace NaN with zeros
df0=df.fillna(value=0)

m = GEKKO()
m.time = df0.index.values

meas1 = m.Param(df0['meas1'].values)
meas2 = m.Param(df0['meas2'].values)
meas3 = m.Param(df0['meas3'].values)
meas4 = m.Param(df0['meas4'].values)
meas5 = m.Param(df0['meas5'].values)
meas6 = m.Param(df0['meas6'].values)

#adjustable Parameters
cnf01=m.FV(1.3,lb=0.01,ub=10) 
cns01=m.FV(1.3,lb=0.01,ub=10) 
kf=m.FV(1.3,lb=0.01,ub=10) 
ks=m.FV(1.3,lb=0.01,ub=10) 

cnf02=m.FV(value=cnf01*0.5,lb=cnf01*0.5, ub=cnf01*0.5)
cns02=m.FV(value=cns01*0.5,lb=cns01*0.5, ub=cns01*0.5)
cnf03=m.FV(value=cnf01*0.25,lb=cnf01*0.25, ub=cnf01*0.25)
cns03=m.FV(value=cns01*0.25,lb=cns01*0.25, ub=cns01*0.25)
cnf04=m.FV(value=cnf01,lb=cnf01, ub=cnf01)
cns04=m.FV(value=cns01,lb=cns01, ub=cns01)
cnf05=m.FV(value=cnf01*0.5,lb=cnf01*0.5, ub=cnf01*0.5)
cns05=m.FV(value=cns01*0.5,lb=cns01*0.5, ub=cns01*0.5)
cnf06=m.FV(value=cnf01*0.25,lb=cnf01*0.25, ub=cnf01*0.25)
cns06=m.FV(value=cns01*0.25,lb=cns01*0.25, ub=cns01*0.25)

#Variables
c1 = m.Var(value=mca1[0])
c2 = m.Var(value=mca2[0])
c3 = m.Var(value=mca3[0])
c4 = m.Var(value=mca4[0])
c5 = m.Var(value=mca5[0])
c6 = m.Var(value=mca6[0])
cm1 = m.Param(df0['x1'].values)
cm2 = m.Param(df0['x2'].values)
cm3 = m.Param(df0['x3'].values)
cm4 = m.Param(df0['x4'].values)
cm5 = m.Param(df0['x5'].values)
cm6 = m.Param(df0['x6'].values)

m.Minimize((meas1*(c1-cm1)**2)+(meas2*(c2-cm2)**2)+(meas3*(c3-cm3)**2)+(meas4*(c4-cm4)**2)+(meas5*(c5-cm5)**2)+(meas6*(c6-cm6)**2))

cnf1=m.Var(value=cnf01)
cns1=m.Var(value=cns01)
cnf2=m.Var(value=cnf02)
cns2=m.Var(value=cns02)
cnf3=m.Var(value=cnf03)
cns3=m.Var(value=cns03)
cnf4=m.Var(value=cnf04)
cns4=m.Var(value=cns04)
cnf5=m.Var(value=cnf05)
cns5=m.Var(value=cns05)
cnf6=m.Var(value=cnf06)
cns6=m.Var(value=cns06)

#Equations
t = m.Param(value=m.time)

m.Equation(cnf1.dt()==-kf*c1*cnf1)
m.Equation(cns1.dt()==-ks*c1*cns1)
m.Equation(c1.dt()==cnf1.dt()+cns1.dt())

m.Equation(cnf2.dt()==-kf*c2*cnf2)
m.Equation(cns2.dt()==-ks*c2*cns2)
m.Equation(c2.dt()==cnf2.dt()+cns2.dt())

m.Equation(cnf3.dt()==-kf*c3*cnf3)
m.Equation(cns3.dt()==-ks*c3*cns3)
m.Equation(c3.dt()==cnf3.dt()+cns3.dt())

m.Equation(cnf4.dt()==-kf*c4*cnf4)
m.Equation(cns4.dt()==-ks*c4*cns4)
m.Equation(c4.dt()==cnf4.dt()+cns4.dt())

m.Equation(cnf5.dt()==-kf*c5*cnf5)
m.Equation(cns5.dt()==-ks*c5*cns5)
m.Equation(c5.dt()==cnf5.dt()+cns5.dt())

m.Equation(cnf6.dt()==-kf*c6*cnf6)
m.Equation(cns6.dt()==-ks*c6*cns6)
m.Equation(c6.dt()==cnf6.dt()+cns6.dt())

m.Equation(ks>0)
m.Equation(kf>0)
m.Equation(cnf01>0)
m.Equation(cns01>0)

#Options
m.options.SOLVER = 3 #IPOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation = MHE
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 3 #collocation nodes (2,5)
m.solve(disp=True)

if True:
    kf.STATUS=1
    ks.STATUS=1
    cnf01.STATUS=1
    cns01.STATUS=1
    cnf02.STATUS=1
    cns02.STATUS=1
    cnf03.STATUS=1
    cns03.STATUS=1
    cnf04.STATUS=1
    cns04.STATUS=1
    cnf05.STATUS=1
    cns05.STATUS=1
    cnf06.STATUS=1
    cns06.STATUS=1

print('Final SSE Objective: ' + str(m.options.objfcnval))

print('Solution')
print('cnf01 = ' + str(cnf01.value[0]))
print('cns01 = ' + str(cns01.value[0]))
print('kf = ' + str(kf.value[0]))
print('ks = ' + str(ks.value[0]))
print('cns02 = '+ str(cns02.value[0]))
print('cnf02 = '+ str(cnf02.value[0]))
print('cns03 = '+ str(cns03.value[0]))
print('cnf03 = '+ str(cnf03.value[0]))
print('cns04 = '+ str(cns04.value[0]))
print('cnf04 = '+ str(cnf04.value[0]))
print('cns05 = '+ str(cns05.value[0]))
print('cnf05 = '+ str(cnf05.value[0]))
print('cns06 = '+ str(cns06.value[0]))
print('cnf06 = '+ str(cnf06.value[0]))

plt.figure(1,figsize=(8,5))
plt.plot(m.time,c1.value,'b',label='Predicted c1')
plt.plot(m.time,c2.value,'r',label='Predicted c2')
plt.plot(m.time,c3.value,'g',label='Predicted c3')
plt.plot(m.time,c4.value,'b',label='Predicted c4')
plt.plot(m.time,c5.value,'r',label='Predicted c5')
plt.plot(m.time,c6.value,'g',label='Predicted c6')
plt.plot(tm1,mca1,'bo',label='Meas c1')
plt.plot(tm2,mca2,'ro',label='Meas c2')
plt.plot(tm3,mca3,'go',label='Meas c3')
plt.plot(tm4,mca4,'bo',label='Meas c4')
plt.plot(tm5,mca5,'ro',label='Meas c5')
plt.plot(tm6,mca6,'go',label='Meas c6')
plt.xlabel('time (h)')
plt.ylabel('Concentration (mg/L)')
plt.legend(loc='best')

Tags： dftimevaluevardtpltvaluesprint
1条回答
网友
1楼 · 发布于 2024-05-23 19:09:44

STATUS=1需要在m.solve()之前，以使kf和ks可由解算器调整

if True:
    kf.STATUS=1
    ks.STATUS=1

不需要额外的不等式约束，因为它们已经在其他地方设置了下限

m.Equation(ks>0)
m.Equation(kf>0)
m.Equation(cnf01>0)
m.Equation(cns01>0)

初始条件通过包括fixed_initial=False来计算

cnf1=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns1=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf2=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns2=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf3=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns3=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf4=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns4=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf5=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns5=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf6=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns6=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)

缩短模型的一种方法是使用数组

cnf1,cnf2,cnf3,cnf4,cnf5,cnf6 = m.Array(m.Var,6,\
            value=1.3,lb=0.0,ub=10,fixed_initial=False
cns1,cns2,cns3,cns4,cns5,cns6 = m.Array(m.Var,6,\
            value=1.3,lb=0.0,ub=10,fixed_initial=False

IPOPT解算器在250次迭代后无法收敛。切换到APOPT解算器将提供成功的解决方案

m.options.SOLVER=1

solution

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt 
import math as math
import pandas as pd

tm1 = [0, 0.0667,0.5,1,4]
mca1 = [5.68, 3.48, 3.24, 3.36, 2.96]

tm2 = [0, 0.08333,0.5,1,4.25]
mca2 = [5.68, 4.20, 4.04, 4.00, 3.76]

tm3 = [0,0.08333,0.5,1,4.33]
mca3 = [5.68, 4.64, 4.52, 4.56, 4.24]

tm4 = [0,0.08333,0.5,1,4.0833]
mca4 =[18.90,15.4,14.3,15.10,13.50]

tm5 = [0,0.08333,0.5,1,4.5]
mca5 =[18.90, 15.5, 16.30, 16, 14.70]

tm6 = [0,0.08333,0.5,1,4.6667]
mca6 = [18.90, 15.8, 11.70,15.5,12]

df1=pd.DataFrame({'time':tm1,'x1':mca1})
df2=pd.DataFrame({'time':tm2,'x2':mca2})
df3=pd.DataFrame({'time':tm3,'x3':mca3})
df4=pd.DataFrame({'time':tm4,'x4':mca4})
df5=pd.DataFrame({'time':tm5,'x5':mca5})
df6=pd.DataFrame({'time':tm6,'x6':mca6})
df1.set_index('time',inplace=True)
df2.set_index('time',inplace=True)
df3.set_index('time',inplace=True)
df4.set_index('time',inplace=True)
df5.set_index('time',inplace=True)
df6.set_index('time',inplace=True)
#simulation time points
dfx = pd.DataFrame({'time':np.linspace(0,5,101)})
dfx.set_index('time',inplace=True)
#merge dataframes
dfxx=dfx.join(df1,how='outer')
dfxxx=dfxx.join(df2,how='outer')
dfxxxx=dfxxx.join(df3,how='outer')
dfxxxxx=dfxxxx.join(df4,how='outer')
dfxxxxxx=dfxxxxx.join(df5,how='outer')
df=dfxxxxxx.join(df6,how='outer')
# get True (1) or False (0) for measurement
df['meas1']=(df['x1'].values==df['x1'].values).astype(int)
df['meas2']=(df['x2'].values==df['x2'].values).astype(int)
df['meas3']=(df['x3'].values==df['x3'].values).astype(int)
df['meas4']=(df['x4'].values==df['x4'].values).astype(int)
df['meas5']=(df['x5'].values==df['x5'].values).astype(int)
df['meas6']=(df['x6'].values==df['x6'].values).astype(int)
#replace NaN with zeros
df0=df.fillna(value=0)

m = GEKKO()
m.time = df0.index.values

meas1 = m.Param(df0['meas1'].values)
meas2 = m.Param(df0['meas2'].values)
meas3 = m.Param(df0['meas3'].values)
meas4 = m.Param(df0['meas4'].values)
meas5 = m.Param(df0['meas5'].values)
meas6 = m.Param(df0['meas6'].values)

#adjustable Parameters
kf=m.FV(1.3,lb=0.01,ub=10) 
ks=m.FV(1.3,lb=0.01,ub=10) 

#Variables
c1 = m.Var(value=mca1[0])
c2 = m.Var(value=mca2[0])
c3 = m.Var(value=mca3[0])
c4 = m.Var(value=mca4[0])
c5 = m.Var(value=mca5[0])
c6 = m.Var(value=mca6[0])
cm1 = m.Param(df0['x1'].values)
cm2 = m.Param(df0['x2'].values)
cm3 = m.Param(df0['x3'].values)
cm4 = m.Param(df0['x4'].values)
cm5 = m.Param(df0['x5'].values)
cm6 = m.Param(df0['x6'].values)

m.Minimize((meas1*(c1-cm1)**2)+(meas2*(c2-cm2)**2)\
          +(meas3*(c3-cm3)**2)+(meas4*(c4-cm4)**2)\
          +(meas5*(c5-cm5)**2)+(meas6*(c6-cm6)**2))

cnf1,cnf2,cnf3,cnf4,cnf5,cnf6 = m.Array(m.Var,6,\
            value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns1,cns2,cns3,cns4,cns5,cns6 = m.Array(m.Var,6,\
            value=1.3,lb=0.0,ub=10,fixed_initial=False)

#Equations
t = m.Param(value=m.time)

m.Equation(cnf1.dt()==-kf*c1*cnf1)
m.Equation(cns1.dt()==-ks*c1*cns1)
m.Equation(c1.dt()==cnf1.dt()+cns1.dt())

m.Equation(cnf2.dt()==-kf*c2*cnf2)
m.Equation(cns2.dt()==-ks*c2*cns2)
m.Equation(c2.dt()==cnf2.dt()+cns2.dt())

m.Equation(cnf3.dt()==-kf*c3*cnf3)
m.Equation(cns3.dt()==-ks*c3*cns3)
m.Equation(c3.dt()==cnf3.dt()+cns3.dt())

m.Equation(cnf4.dt()==-kf*c4*cnf4)
m.Equation(cns4.dt()==-ks*c4*cns4)
m.Equation(c4.dt()==cnf4.dt()+cns4.dt())

m.Equation(cnf5.dt()==-kf*c5*cnf5)
m.Equation(cns5.dt()==-ks*c5*cns5)
m.Equation(c5.dt()==cnf5.dt()+cns5.dt())

m.Equation(cnf6.dt()==-kf*c6*cnf6)
m.Equation(cns6.dt()==-ks*c6*cns6)
m.Equation(c6.dt()==cnf6.dt()+cns6.dt())

if True:
    kf.STATUS=1
    ks.STATUS=1

#Options
m.options.SOLVER = 1 #IPOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation = MHE
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 3 #collocation nodes (2,5)
m.solve(disp=True)

print('Final SSE Objective: ' + str(m.options.objfcnval))

print('Solution')
print('kf = ' + str(kf.value[0]))
print('ks = ' + str(ks.value[0]))

if True:
    print('cnf01 = ' + str(cnf1.value[0]))
    print('cns01 = ' + str(cns1.value[0]))
    print('cns02 = '+ str(cns2.value[0]))
    print('cnf02 = '+ str(cnf2.value[0]))
    print('cns03 = '+ str(cns3.value[0]))
    print('cnf03 = '+ str(cnf3.value[0]))
    print('cns04 = '+ str(cns4.value[0]))
    print('cnf04 = '+ str(cnf4.value[0]))
    print('cns05 = '+ str(cns5.value[0]))
    print('cnf05 = '+ str(cnf5.value[0]))
    print('cns06 = '+ str(cns6.value[0]))
    print('cnf06 = '+ str(cnf6.value[0]))

plt.figure(1,figsize=(8,5))
plt.plot(m.time,c1.value,'b',label='Predicted c1')
plt.plot(m.time,c2.value,'r',label='Predicted c2')
plt.plot(m.time,c3.value,'g',label='Predicted c3')
plt.plot(m.time,c4.value,'b',label='Predicted c4')
plt.plot(m.time,c5.value,'r',label='Predicted c5')
plt.plot(m.time,c6.value,'g',label='Predicted c6')
plt.plot(tm1,mca1,'bo',label='Meas c1')
plt.plot(tm2,mca2,'ro',label='Meas c2')
plt.plot(tm3,mca3,'go',label='Meas c3')
plt.plot(tm4,mca4,'bo',label='Meas c4')
plt.plot(tm5,mca5,'ro',label='Meas c5')
plt.plot(tm6,mca6,'go',label='Meas c6')
plt.xlabel('time (h)')
plt.ylabel('Concentration (mg/L)')
plt.legend(loc='best')
plt.show()

                          -
 Solver         :  APOPT (v1.0)
 Solution time  :    20.3781000000017      sec
 Objective      :    33.5882065490541     
 Successful solution
                          -
 
Final SSE Objective: 33.588206549
Solution
kf = 0.01
ks = 1.0397291201
cnf01 = 0.0
cns01 = 2.7570924931
cns02 = 1.9055347775
cnf02 = 0.0
cns03 = 1.2943612345
cnf03 = 0.58537442202
cns04 = 3.9482253816
cnf04 = 3.0425925867
cns05 = 2.8833138483
cnf05 = 2.3309239496
cns06 = 4.6060598019
cnf06 = 4.5472709683

相关问题 更多 >

    编程相关推荐