问题描述
解决问题 How to apply m.connection from Gekko using arrays?
我尝试使用新策略:
我开始测试是否可以在估计之前进行模拟,因为我发现了一些用作参数的值。为简化起见,我仅模拟一个数据,并且通过假设参数 kdoc1 和 kdoc2、doc1(t=0) 和 doc2(t=0) 为零,我假设氨是唯一影响变量 C2 (NH2Cl) 衰减的反应性试剂,所以 r15 和 r16 也为零。这里的结果是有道理的,因为衰变只是由于氨。为了获得更高的衰减和更好的数据近似值(红色方块),r15 和 r16 不应为零,这也表明有机物在 C2 (NH2Cl) 的衰减中的贡献。
这是图表结果:
然后我的想法是进行模拟,为参数添加值(这些值经过验证)。即使 DOF=0 我也有错误 1. 是否有可能我需要使用 m.connection 来通知 doc10 is doc1 at t=0 即使在模拟模式下?或者需要将变量标识为SV。还是FV。?我尝试申请,但没有成功。提前致谢。
这是我的代码,假设参数 kdoc1、kdoc2、doc10 (doc1(t=0)) 和 doc20 (doc2(t=0)) 的值:
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)
df['meas'] = (df['xmgl'].values==df['xmgl'].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)
nh3 = m.Var(value=1.37e-3)
nh2cl = m.Var(value=0)
nhcl2 = m.Var(value=0)
h = m.Var(value=6.37e-8)
oh = m.Var(value=1e-14/h)
I = m.Var(value=1e-15)
ocl = m.Var(value=0)
nh4 = m.Var(value=0)
h2co3 = m.Var(value=5.39e-4)
hco3 = m.Var(value=3.46e-3)
co32 = m.Var(value=2.16e-6)
alk = m.Var(value=3.46e-3)
###adjustable Parameters
doc10 = 0 #doc1(t=0)
doc20 = 0 #doc2(t=0)
kdoc1 = 0
kdoc2 = 0
#doc10 = 3.1684e-5
#doc20 = 3.1308e-5
#kdoc1 = 1.2334e6
#kdoc2 = 3.8809e9
doc1 = m.Var(value=doc10)
doc2 = m.Var(value=doc20)
####Rate constants
k1 = 1.5e10
k2 = 7.6e-2
k3 = 1e6
k4 = 2.3e-3
k5 = 2.5e7*h+4e4*h2co3+800*hco3
k6 = 2.2e8
k7 = 4e5
k8 = 1e8
k9 = 3e7
k10 = 55
k11=3.16e-8*1e10
k12=1e10
k13=5.01e-10*1e10
k14=1e10
r1 = k1 * hocl * nh3
r2 = k2 * nh2cl
r3 = k3 * hocl * nh2cl
r4 = k4 * nhcl2
r5 = k5 * nh2cl * nh2cl
r6 = k6 * nhcl2 * nh3* h
r7 = k7 * nhcl2 * oh
r8 = k8 * I * nhcl2
r9 = k9 * I * nh2cl
r10 = k10 * nh2cl * nhcl2
r11=k11*hocl
r12=k12*h*ocl
r13=k13*nh4
r14=k14*h*nh3
r15 = kdoc1*doc1*nh2cl
r16 = kdoc2*doc2*hocl
t = m.Param(value=m.time)
#m.Connection(doc1,doc10,pos1=1,pos2=1,node1=1,node2=1)
#m.Connection(doc2,doc20,node2=1)
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)
m.options.IMODE=4 #dynamic simultaneous / Simulation
#m.options.NODES=2
m.solve()
#m.solve(disp=False)
#apm_get(server,app,'infeasibilities.txt')
#m.open_folder()
###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,2)
plt.plot(m.time,nh3.value,label ='NH3')
plt.legend()
plt.subplot(4,3)
plt.plot(m.time,label ='NHCl2')
plt.legend()
plt.subplot(4,4)
plt.plot(m.time,h.value,label ='H+')
plt.legend()
plt.subplot(4,5)
plt.plot(m.time,oh.value,label ='OH-')
plt.legend()
plt.subplot(4,6)
plt.plot(m.time,I.value,label ='I')
plt.legend()
plt.subplot(4,7)
plt.plot(m.time,ocl.value,label ='OCl-')
plt.legend()
plt.subplot(4,8)
plt.plot(m.time,nh4.value,label ='NH4+')
plt.legend()
plt.subplot(4,9)
plt.plot(m.time,h2co3.value,label ='H2CO3')
plt.legend()
plt.xlabel('time (h)')
plt.subplot(4,10)
plt.plot(m.time,hco3.value,label ='HCO3-')
plt.legend()
plt.xlabel('time (h)')
plt.subplot(4,11)
plt.plot(m.time,co32.value,label ='CO32-')
plt.legend()
plt.xlabel('time (h)')
plt.subplot(4,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)'])
解决方法
要回答如何设置dynamic的初始条件的问题,value
是在模拟模式下指定初始条件的正确方法,例如x=m.Var(value=0)
。
doc1 = m.Var(value=doc10)
如果计算了初始条件,则将 fixed_initial=False
设置为 IMODE=5
以允许自由度。
m.options.IMODE=5
doc1 = m.Var(value=doc10,fixed_initial=False)
doc10
值给出了初始猜测值,然后优化器会更改该值以最小化目标函数,例如拟合数据的平方误差。
使用 IMODE=7
产生与 IMODE=4
等效的结果,但按顺序而不是同时求解动态方程。通常更容易解决。