微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦！

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

分类：编程问答

如何解决使用 GEKKO二阶微分方程估计参数

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

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,4.25]
mca2 = [5.68,4.20,4.04,4.00,3.76]

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

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

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

tm6 = [0,4.6667]
mca6 = [18.90,15.8,11.70,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,ub=10) 
kf=m.FV(1.3,ub=10) 
ks=m.FV(1.3,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,ub=cnf01*0.5)
cns05=m.FV(value=cns01*0.5,ub=cns01*0.5)
cnf06=m.FV(value=cnf01*0.25,ub=cnf01*0.25)
cns06=m.FV(value=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,label='Predicted c4')
plt.plot(m.time,c5.value,label='Predicted c5')
plt.plot(m.time,c6.value,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,label='Meas c4')
plt.plot(tm5,mca5,label='Meas c5')
plt.plot(tm6,mca6,label='Meas c6')
plt.xlabel('time (h)')
plt.ylabel('Concentration (mg/L)')
plt.legend(loc='best')

解决方法

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,fixed_initial=False)
cnf2=m.Var(value=1.3,fixed_initial=False)
cns2=m.Var(value=1.3,fixed_initial=False)
cnf3=m.Var(value=1.3,fixed_initial=False)
cns3=m.Var(value=1.3,fixed_initial=False)
cnf4=m.Var(value=1.3,fixed_initial=False)
cns4=m.Var(value=1.3,fixed_initial=False)
cnf5=m.Var(value=1.3,fixed_initial=False)
cns5=m.Var(value=1.3,fixed_initial=False)
cnf6=m.Var(value=1.3,fixed_initial=False)
cns6=m.Var(value=1.3,fixed_initial=False)

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

cnf1,cnf2,cnf3,cnf4,cnf5,cnf6 = m.Array(m.Var,6,\
            value=1.3,fixed_initial=False
cns1,cns2,cns3,cns4,cns5,cns6 = m.Array(m.Var,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,4.25]
mca2 = [5.68,4.20,4.04,4.00,3.76]

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

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

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

tm6 = [0,4.6667]
mca6 = [18.90,15.8,11.70,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,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,fixed_initial=False)
cns1,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,label='Predicted c4')
plt.plot(m.time,c5.value,label='Predicted c5')
plt.plot(m.time,c6.value,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,label='Meas c4')
plt.plot(tm5,mca5,label='Meas c5')
plt.plot(tm6,mca6,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

版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。

苹果市值2025年有望达4万亿美元
pythonJavaScriptjavaHTMLPHPreactjsC#AndroidCSSNode.jssqlrpython-3.xMysqLjQueryc++pandasFlutterangularIOSdjangolinuxswifttypescript路由器JSON路由器设置无线路由器h3c华三华三路由器设置华三路由器电脑软件教程arraysdocker软件图文教程Cvue.jslaravelspring-boot