利用gekko关于模型结构的建议估算速率系数和初始浓度

Question

现在我使用的是简化代码，在这里我得到了一个模拟结果（没有数组，只针对一个数据集）： Do I need to use Gekko m.connection in simulation mode to solve my code?

我正试图将其更改为评估/优化模型，与下面的模型相对应。极小化函数是否和m连接一样好？有什么不正确的建议吗？先谢谢你

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

###measured data1 - nh2cl
data_mgl = [0,3.48,3.24,3.36,2.96,1.96]#5.68  
data =  [0,4.91e-5,4.57e-5,4.74e-5,4.17e-5,2.76e-5] #8.01e-5
t_data = [0,0.08333,0.5,1,4,22.6167]

df1 = pd.DataFrame({'time':t_data,'x':data,'xmgl':data_mgl})
df1.set_index('time', inplace=True)
# simulation time points
df2 = pd.DataFrame({'time':np.linspace(0,25,200)}) #120
df2.set_index('time', inplace=True)
# merge dataframes
df = df2.join(df1,how='outer')
# get True (1) or False (0) for measurement
df['meas'] = (df['x'].values==df['x'].values).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)

###Estimator Model

m = GEKKO(remote=False)

m.time = df0.index.values

#species concentration (M=mol/l)

hocl = m.Var(value=8.01e-5,lb=0)
nh3 = m.Var(value=1.37e-3,lb=0)
nh2cl = m.Var(value=0,ub=5.5e-5,lb=0)
nhcl2 = m.Var(value=0,lb=0)
h = m.Var(value=6.37e-8,lb=0)
oh = m.Var(value=1e-14/h,lb=0)
I = m.Var(value=1e-15,lb=0) 
ocl = m.Var(value=0,lb=0)
nh4 = m.Var(value=0,lb=0)
h2co3 = m.Var(value=5.39e-4,lb=0) 
hco3 = m.Var(value=3.46e-3,lb=0)
co32 = m.Var(value=2.16e-6,lb=0) 
alk = m.Var(value=3.46e-3,lb=0) 

####Rate constants
k1 = m.Const(1.5e10)
k2 = m.Const(7.6e-2)
k3 = m.Const(1e6)
k4 = m.Const(2.3e-3)
k6 = m.Const(2.2e8)
k7 = m.Const(4e5)
k8 = m.Const(1e8)
k9 = m.Const(3e7)
k10 = m.Const(55)

k11=m.Const(3.16e-8*1e10)
k12=m.Const(1e10)
k13=m.Const(5.01e-10*1e10)
k14=m.Const(1e10)

k5 =m.Var(2.5e7*h+4e4*h2co3+800*hco3)

r1=m.Var()
r2=m.Var()
r3=m.Var()
r4=m.Var()
r5=m.Var()
r6=m.Var()
r7=m.Var()
r8=m.Var()
r9=m.Var()
r10=m.Var()
r11=m.Var()
r12=m.Var()
r13=m.Var()
r14=m.Var()
r15=m.Var()
r16=m.Var()

###Adjustable Parameters

kdoc1 = m.FV(1.2334e6,lb=1,ub=1e10) 
kdoc2 = m.FV(3.8809e9,lb=1,ub=1e10)
DOC10 = m.FV(3.1684e-5,lb=1e-6,ub=1)
DOC20 = m.FV(3.1308e-5,lb=1e-6,ub=1)

#doc10 = 3.1684e-5
#doc20 = 3.1308e-5
#kdoc1 = 1.2334e6 
#kdoc2 = 3.8809e9

doc1 = m.SV(value=DOC10)
doc2 = m.SV(value=DOC20)

m.Connection(doc1,DOC10,pos1=1,pos2=1,node1=1,node2=1)
m.Connection(doc2,DOC20,pos1=1,pos2=1,node1=1,node2=1)

m.Equation(k5 == 2.5e7*h+4e4*h2co3+800*hco3)

m.Equation(r1 == k1 * hocl * nh3)
m.Equation(r2 == k2 * nh2cl)
m.Equation(r3 == k3 * hocl * nh2cl)
m.Equation(r4 == k4 * nhcl2)
m.Equation(r5 == k5 * nh2cl * nh2cl)
m.Equation(r6 == k6 * nhcl2 * nh3* h)
m.Equation(r7 == k7 * nhcl2 * oh)
m.Equation(r8 == k8 * I * nhcl2)
m.Equation(r9 == k9 * I * nh2cl)
m.Equation(r10 == k10 * nh2cl * nhcl2)

m.Equation(r11 == k11*hocl)
m.Equation(r12 == k12*h*ocl)
m.Equation(r13 == k13*nh4)
m.Equation(r14 == k14*h*nh3)

m.Equation(r15 == kdoc1*doc1*nh2cl)
m.Equation(r16 == kdoc2*doc2*hocl)

m.Equation(co32 == (5.01e-11*hco3)/h)
m.Equation(h2co3 == (hco3*h)/5.01e-7)
m.Equation(hco3 == alk - 2*co32 - oh + h) 
m.Equation(oh == 1e-14/h)

m.Equation(hocl.dt()== -r1 + r2 - r3 + r4 + r8 - r11 + r12 - r16)
m.Equation(nh3.dt()== -r1 + r2 + r5 - r6 + r13 - r14)
m.Equation(nh2cl.dt()== r1 - r2 - r3 + r4 - r5 + r6 - r9 - r10 - r15)
m.Equation(nhcl2.dt()== r3 - r4 + r5 - r6 - r7 - r8 - r10)
m.Equation(h.dt()== 0)
m.Equation(I.dt()== r7 - r8 - r9)
m.Equation(alk.dt()== 0)
m.Equation(ocl.dt()==r11-r12)
m.Equation(nh4.dt()==-r13+r14)

m.Equation(doc1.dt()==-r15)
m.Equation(doc2.dt()==-r16)

meas1 = m.Param(df0['meas'].values)
cm1 = m.Param(df0['x'].values)

#m.Minimize(meas1*((nh2cl-df0['x'].values))**2)
m.Minimize(meas1*((nh2cl-cm1))**2)

m.options.SOLVER = 1
m.options.IMODE = 8 
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 3
m.options.RTOL = 1E-3

if True:
    kdoc1.STATUS=1
    kdoc2.STATUS=1
    DOC10.STATUS=1
    DOC20.STATUS=1

m.solve()
#m.solve(disp=False)
m.open_folder()

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

print('Solution')
print('kdoc1 = ' + str(kdoc1))
print('kdoc2 = ' + str(kdoc2))
print('DOC10 = ' + str(DOC10))
print('DOC20 = ' + str(DOC20))

###Graphics

plt.xlabel('time (h)')
plt.ylabel('Concentration (mg/L)')
plt.legend(loc='best')

plt.figure(2)
#plt.plot(m.time,cC2.value,label ='NH2Cl')
#plt.plot(m.time,cC1.value,label ='NH3')
#plt.plot(m.time,C0.value,label ='HOCl')
#plt.plot(m.time,C1.value,label ='NH3')
plt.plot(m.time,nh2cl.value,label ='NH2Cl')
plt.plot(m.time,nhcl2.value,label ='NHCl2')
#plt.plot(m.time,C4.value,label ='H+')
#plt.plot(m.time,C5.value,label ='OH-')
#plt.plot(m.time,C6.value,label ='I')
#plt.plot(m.time,C7.value,label ='OCl-')
#plt.plot(m.time,C8.value,label ='NH4+')
#plt.plot(m.time,C9.value,label ='H2CO3')
#plt.plot(m.time,C10.value,label ='HCO3-')
#plt.plot(m.time,C11.value,label ='CO2-3')
#plt.plot(m.time,C12.value,label ='alk')
#plt.plot(m.time,TotalCl.value,label ='TotalCl')
plt.plot(m.time,df['x'].values,'rs',label='Meas')

plt.xlabel('time (h)')
plt.ylabel('Concentration (mol/L)')
plt.legend(loc='best')
plt.grid()
plt.xlim(-0.05, 25)

plt.figure(3,figsize=(12,8))
plt.subplot(4,3,1)
plt.plot(m.time,hocl.value,label ='HOCl')
plt.legend()

plt.subplot(4,3,2)
plt.plot(m.time,nh3.value,label ='NH3')
plt.legend()

plt.subplot(4,3,3)
plt.plot(m.time,nhcl2.value,label ='NHCl2')
plt.legend()

plt.subplot(4,3,4)
plt.plot(m.time,h.value,label ='H+')
plt.legend()

plt.subplot(4,3,5)
plt.plot(m.time,oh.value,label ='OH-')
plt.legend()

plt.subplot(4,3,6)
plt.plot(m.time,I.value,label ='I')
plt.legend()

plt.subplot(4,3,7)
plt.plot(m.time,ocl.value,label ='OCl-')
plt.legend()

plt.subplot(4,3,8)
plt.plot(m.time,nh4.value,label ='NH4+')
plt.legend()

plt.subplot(4,3,9)
plt.plot(m.time,h2co3.value,label ='H2CO3')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,10)
plt.plot(m.time,hco3.value,label ='HCO3-')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,11)
plt.plot(m.time,co32.value,label ='CO32-')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,12)
plt.plot(m.time,alk.value,label ='alk')
plt.legend()
plt.xlabel('time (h)')

plt.show()

datasa1={'Time (h)':m.time, 'HOCl (M)':hocl.value, 'NH3 (M)':nh3.value, 'NH2Cl (M)':nh2cl.value, 'NHCl2 (M)':nhcl2.value,'H+ (M)':h.value,'OH- (M)':oh.value, 'I (M)':I.value,'OCl- (M)':ocl,'NH4+ (M)':nh4,'H2CO3 (M)':h2co3,'HCO3 (M)':hco3,'CO3 (M)':co32,'Alk (M)':alk}    

dftestsa1=pd.DataFrame(datasa1, columns=['Time (h)', 'HOCl (M)','NH3 (M)','NH2Cl (M)',\
                               'NHCl2 (M)','H+ (M)','OH- (M)',\
                               'I (M)','OCl- (M)','NH4+ (M)','H2CO3 (M)','HCO3 (M)','CO3 (M)',\
                               'Alk (M)'])

Answer 1

估计问题有两个自由度，但要计算的值应该有4个自由度

一些建议：

• 可以使用数组声明变量，以简化声明

r1,r2,r3,r4,r5,r6,r7,r8,\
 r9,r10,r11,r12,r13,r14,r15,r16=m.Array(m.Var,16)

• 更好的方法是将速率常数声明为中间产物

r1 = m.Intermediate(k1 * hocl * nh3)
r2 = m.Intermediate(k2 * nh2cl)
r3 = m.Intermediate(k3 * hocl * nh2cl)
r4 = m.Intermediate(k4 * nhcl2)
r5 = m.Intermediate(k5 * nh2cl * nh2cl)
r6 = m.Intermediate(k6 * nhcl2 * nh3* h)
r7 = m.Intermediate(k7 * nhcl2 * oh)
r8 = m.Intermediate(k8 * I * nhcl2)
r9 = m.Intermediate(k9 * I * nh2cl)
r10 = m.Intermediate(k10 * nh2cl * nhcl2)

r11 = m.Intermediate(k11*hocl)
r12 = m.Intermediate(k12*h*ocl)
r13 = m.Intermediate(k13*nh4)
r14 = m.Intermediate(k14*h*nh3)

r15 = m.Intermediate(kdoc1*doc1*nh2cl)
r16 = m.Intermediate(kdoc2*doc2*hocl)

• 1e-3处的残余公差相当高。我建议将默认值设置为1e-6，特别是对于h和oh这样的小浓度

m.options.RTOL = 1E-6

• 同时模式IMODE=5的效果比顺序估计模式IMODE = 8好得多

• 到DOC10和DOC20的连接是计算初始条件的一种方法。一般来说，最好对DOC1/DOC2使用fixed_initial=False，并对变量施加约束以限制初始条件。但是，此方法应该有效，但自由度仅为2（应为4）

• 我试图用IMODE=4设置初始化。这通常通过在优化之前找到初始可行解来帮助解算器

• 请检查以下方程式。它们暗示hco3==alk - 2*5.01e-11*hco3/h - oh + h具有替换和代数重排。这是正确的吗

m.Equation(co32 == (5.01e-11*hco3)/h)
m.Equation(h2co3 == (hco3*h)/5.01e-7)
m.Equation(hco3 == alk - 2*co32 - oh + h) 
m.Equation(oh == 1e-14/h)

下面是脚本的当前版本，该版本尚未运行，但更接近于提供解决方案。请让我们了解上述方程式

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

###measured data1 - nh2cl
data_mgl = [0,3.48,3.24,3.36,2.96,1.96]#5.68  
data =  [0,4.91e-5,4.57e-5,4.74e-5,4.17e-5,2.76e-5] #8.01e-5
t_data = [0,0.08333,0.5,1,4,22.6167]

df1 = pd.DataFrame({'time':t_data,'x':data,'xmgl':data_mgl})
df1.set_index('time', inplace=True)
# simulation time points
df2 = pd.DataFrame({'time':np.linspace(0,25,200)}) #120
df2.set_index('time', inplace=True)
# merge dataframes
df = df2.join(df1,how='outer')
# get True (1) or False (0) for measurement
df['meas'] = (df['x'].values==df['x'].values).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)

###Estimator Model

m = GEKKO(remote=True)

m.time = df0.index.values

#species concentration (M=mol/l)

hocl = m.Var(value=8.01e-5,lb=0)
nh3 = m.Var(value=1.37e-3,lb=0)
nh2cl = m.Var(value=0,ub=5.5e-5,lb=0)
nhcl2 = m.Var(value=0,lb=0)

# equation for h is not needed
h = m.Param(value=6.37e-8)# ,lb=0)
#m.Equation(h.dt()== 0)

# use Intemediate variable for OH concentration
#oh = m.Var(value=1e-14/6.37e-8,lb=0)
#m.Equation(oh * h == 1e-14)
oh = m.Intermediate(1e-14/h)

I = m.Var(value=1e-15,lb=0) 
ocl = m.Var(value=0,lb=0)
nh4 = m.Var(value=0,lb=0)
#h2co3 = m.Var(value=5.39e-4,lb=0) 
#co32 = m.Var(value=2.16e-6,lb=0) 
alk = m.Var(value=3.46e-3,lb=0)

hco3  = m.Var(3.46e-3,lb=0)
m.Equation(hco3==alk - 2*5.01e-11*hco3/h - oh + h) 
co32  = m.Intermediate(5.01e-11*hco3/h)
h2co3 = m.Intermediate((hco3*h)/5.01e-7)

####Rate constants
k1 = m.Const(1.5e10)
k2 = m.Const(7.6e-2)
k3 = m.Const(1e6)
k4 = m.Const(2.3e-3)
k6 = m.Const(2.2e8)
k7 = m.Const(4e5)
k8 = m.Const(1e8)
k9 = m.Const(3e7)
k10 = m.Const(55)
k11=m.Const(3.16e-8*1e10)
k12=m.Const(1e10)
k13=m.Const(5.01e-10*1e10)
k14=m.Const(1e10)

k5 =m.Var(2.5e7*h+4e4*h2co3+800*hco3)

#r1,r2,r3,r4,r5,r6,r7,r8,\
# r9,r10,r11,r12,r13,r14,r15,r16=m.Array(m.Var,16)

###Adjustable Parameters

kdoc1 = m.FV(1.2334e6,lb=1,ub=1e10) 
kdoc2 = m.FV(3.8809e9,lb=1,ub=1e10)
DOC10 = m.FV(3.1684e-5,lb=1e-6,ub=1)
DOC20 = m.FV(3.1308e-5,lb=1e-6,ub=1)

#doc10 = 3.1684e-5
#doc20 = 3.1308e-5
#kdoc1 = 1.2334e6 
#kdoc2 = 3.8809e9

doc1 = m.SV(value=DOC10)
doc2 = m.SV(value=DOC20)

m.Equation(k5 == 2.5e7*h+4e4*h2co3+800*hco3)

r1 = m.Intermediate(k1 * hocl * nh3)
r2 = m.Intermediate(k2 * nh2cl)
r3 = m.Intermediate(k3 * hocl * nh2cl)
r4 = m.Intermediate(k4 * nhcl2)
r5 = m.Intermediate(k5 * nh2cl * nh2cl)
r6 = m.Intermediate(k6 * nhcl2 * nh3* h)
r7 = m.Intermediate(k7 * nhcl2 * oh)
r8 = m.Intermediate(k8 * I * nhcl2)
r9 = m.Intermediate(k9 * I * nh2cl)
r10 = m.Intermediate(k10 * nh2cl * nhcl2)

r11 = m.Intermediate(k11*hocl)
r12 = m.Intermediate(k12*h*ocl)
r13 = m.Intermediate(k13*nh4)
r14 = m.Intermediate(k14*h*nh3)

r15 = m.Intermediate(kdoc1*doc1*nh2cl)
r16 = m.Intermediate(kdoc2*doc2*hocl)

m.Equation(hocl.dt()== -r1 + r2 - r3 + r4 + r8 - r11 + r12 - r16)
m.Equation(nh3.dt()== -r1 + r2 + r5 - r6 + r13 - r14)
m.Equation(nh2cl.dt()== r1 - r2 - r3 + r4 - r5 + r6 - r9 - r10 - r15)
m.Equation(nhcl2.dt()== r3 - r4 + r5 - r6 - r7 - r8 - r10)
m.Equation(I.dt()== r7 - r8 - r9)
m.Equation(alk.dt()== 0)
m.Equation(ocl.dt()==r11-r12)
m.Equation(nh4.dt()==-r13+r14)

m.Equation(doc1.dt()==-r15)
m.Equation(doc2.dt()==-r16)

meas1 = m.Param(df0['meas'].values)
cm1 = m.Param(df0['x'].values)

#m.Minimize(meas1*((nh2cl-df0['x'].values))**2)
m.Minimize(meas1*((nh2cl-cm1))**2)

m.options.SOLVER = 1
m.options.EV_TYPE = 2 #squared error
m.options.NODES = 3
m.options.RTOL = 1E-6

if True:
    kdoc1.STATUS=1
    kdoc2.STATUS=1
    DOC10.STATUS=1
    DOC20.STATUS=1

m.options.IMODE = 4
m.options.RTOL = 1e-5
m.solve()

m.Connection(doc1,DOC10,pos1=1,pos2=1,node1=1,node2=1)
m.Connection(doc2,DOC20,pos1=1,pos2=1,node1=1,node2=1)

m.options.IMODE = 5
m.options.COLDSTART=0
m.options.TIME_SHIFT=0
m.solve()
#m.solve(disp=False)
m.open_folder()

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

print('Solution')
print('kdoc1 = ' + str(kdoc1))
print('kdoc2 = ' + str(kdoc2))
print('DOC10 = ' + str(DOC10))
print('DOC20 = ' + str(DOC20))

###Graphics

plt.xlabel('time (h)')
plt.ylabel('Concentration (mg/L)')
plt.legend(loc='best')

plt.figure(2)
#plt.plot(m.time,cC2.value,label ='NH2Cl')
#plt.plot(m.time,cC1.value,label ='NH3')
#plt.plot(m.time,C0.value,label ='HOCl')
#plt.plot(m.time,C1.value,label ='NH3')
plt.plot(m.time,nh2cl.value,label ='NH2Cl')
plt.plot(m.time,nhcl2.value,label ='NHCl2')
#plt.plot(m.time,C4.value,label ='H+')
#plt.plot(m.time,C5.value,label ='OH-')
#plt.plot(m.time,C6.value,label ='I')
#plt.plot(m.time,C7.value,label ='OCl-')
#plt.plot(m.time,C8.value,label ='NH4+')
#plt.plot(m.time,C9.value,label ='H2CO3')
#plt.plot(m.time,C10.value,label ='HCO3-')
#plt.plot(m.time,C11.value,label ='CO2-3')
#plt.plot(m.time,C12.value,label ='alk')
#plt.plot(m.time,TotalCl.value,label ='TotalCl')
plt.plot(m.time,df['x'].values,'rs',label='Meas')

plt.xlabel('time (h)')
plt.ylabel('Concentration (mol/L)')
plt.legend(loc='best')
plt.grid()
plt.xlim(-0.05, 25)

plt.figure(3,figsize=(12,8))
plt.subplot(4,3,1)
plt.plot(m.time,hocl.value,label ='HOCl')
plt.legend()

plt.subplot(4,3,2)
plt.plot(m.time,nh3.value,label ='NH3')
plt.legend()

plt.subplot(4,3,3)
plt.plot(m.time,nhcl2.value,label ='NHCl2')
plt.legend()

plt.subplot(4,3,4)
plt.plot(m.time,h.value,label ='H+')
plt.legend()

plt.subplot(4,3,5)
plt.plot(m.time,oh.value,label ='OH-')
plt.legend()

plt.subplot(4,3,6)
plt.plot(m.time,I.value,label ='I')
plt.legend()

plt.subplot(4,3,7)
plt.plot(m.time,ocl.value,label ='OCl-')
plt.legend()

plt.subplot(4,3,8)
plt.plot(m.time,nh4.value,label ='NH4+')
plt.legend()

plt.subplot(4,3,9)
plt.plot(m.time,h2co3.value,label ='H2CO3')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,10)
plt.plot(m.time,hco3.value,label ='HCO3-')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,11)
plt.plot(m.time,co32.value,label ='CO32-')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,12)
plt.plot(m.time,alk.value,label ='alk')
plt.legend()
plt.xlabel('time (h)')

plt.show()

datasa1={'Time (h)':m.time, 'HOCl (M)':hocl.value, 'NH3 (M)':nh3.value, 'NH2Cl (M)':nh2cl.value, 'NHCl2 (M)':nhcl2.value,'H+ (M)':h.value,'OH- (M)':oh.value, 'I (M)':I.value,'OCl- (M)':ocl,'NH4+ (M)':nh4,'H2CO3 (M)':h2co3,'HCO3 (M)':hco3,'CO3 (M)':co32,'Alk (M)':alk}    

dftestsa1=pd.DataFrame(datasa1, columns=['Time (h)', 'HOCl (M)','NH3 (M)','NH2Cl (M)',\
                               'NHCl2 (M)','H+ (M)','OH- (M)',\
                               'I (M)','OCl- (M)','NH4+ (M)','H2CO3 (M)','HCO3 (M)','CO3 (M)',\
                               'Alk (M)'])