python - 使用 Gekko 的模拟模式 - 初始值为零的变量 C2 (NH2Cl) 不起作用 - 为什么？

Question

如果我假设变量 C2(NH2Cl) 的初始值为零，C2 = m.Var(value=0) 或更接近于零，如 1e-10，则我的模型模拟不会运行。如果 C2 = m.Var(value=0.3e-5)，则模型模拟运行。你知道为什么会这样吗？我该如何解决？先感谢您。

这是 C2 = m.Var(value=0.3e-5) 的图形结果：

这是我的代码：

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

#measured data run1 W1-M1
data_mgl=[0,2.44,1.72,2.32,2.16]
data=[0,3.44e-5,2.43e-5,3.27e-5,3.05e-5]
t_data=[0,0.0833333,0.5,1,4]

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,5,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()
#m = GEKKO(remote=False)
m.time = df0.index.values


#species concentration (M=mol/l)
C0 = m.Var(value=3.05e-5)
C1 = m.Var(value=3.57e-4)
C2 = m.Var(value=0)
#C2 = m.Var(value=0.3e-5)
C3 = m.Var(value=1e-15)
C4 = m.Var(value=1.479e-8)
C5 = m.Var(value=6.761e-7)
C6 = m.Var(value=0)
C7 = m.Var(value=0)
C8 = m.Var(value=0)
C9 = m.Var(value=1.2e-4)
C10 = m.Var(value=3.28e-3)
C11 = m.Var(value=8.72e-6)
C12 = m.Var(value=3.3e-3)
doc1 = m.Var(value=0)
doc2 = m.Var(value=0)

#cC1 = m.Var(value=0)
#cC2 = m.Var(value=0)

#Rate constants
k1 = 1.5e10
k2 = 7.6e-2
k3 = 1e6
k4 = 2.3e-3
k5 = 2.5e7*C4+4e4*C9+800*C10
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

kdoc1 = 0   
kdoc2 = 0     

r1 = k1 * C0 * C1
r2 = k2 * C2
r3 = k3 * C0 * C2
r4 = k4 * C3
r5 = k5 * C2 * C2
r6 = k6 * C3 * C1* C4
r7 = k7 * C3 * C5
r8 = k8 * C6 * C3
r9 = k9 * C6 * C2
r10 = k10 * C2 * C3

r11=k11*C0
r12=k12*C4*C7
r13=k13*C8
r14=k14*C4*C1

r15 = kdoc1*doc1*C2
r16 = kdoc2*doc2*C0

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

m.Equation(C0.dt()== -r1 + r2 - r3 + r4 + r8 - r16 - r11 + r12)
######m.Equation(C1.dt()== -r1 + r2 + r5 - r6 + r15)-porque é retirada a r15?!!!!!!
m.Equation(C1.dt()== -r1 + r2 + r5 - r6 + r13 - r14)
m.Equation(C2.dt()== r1 - r2 - r3 + r4 - r5 + r6 - r9 - r10 - r15)
m.Equation(C3.dt()== r3 - r4 + r5 - r6 - r7 - r8 - r10)
m.Equation(C4.dt()== 0)
m.Equation(C5 == 1e-14/C4)
m.Equation(C6.dt()== r7 - r8 - r9)
m.Equation(C7.dt()==r11-r12)
m.Equation(C8.dt()==-r13+r14)
m.Equation(C9 == (C10*C4)/5.01e-7)
m.Equation(C10 == C12 - 2*C11 - C5 + C4) 
m.Equation(C11 == (5.01e-11*C10)/C4)
m.Equation(C12.dt()== 0)

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

m.options.IMODE=4 #dynamic simultaneous / Simulation
m.solve(disp=False)
#apm_get(server,app,'infeasibilities.txt') 
#m.open_folder()

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,C2.value,label ='NH2Cl')
plt.plot(m.time,C3.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, 5)

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

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

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

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

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

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

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

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

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

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

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

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

plt.show()

datasa1={'Time (h)':m.time, 'HOCl (M)':C0.value, 'NH3 (M)':C1.value, 'NH2Cl (M)':C2.value, 'NHCl2 (M)':C3.value,'H+ (M)':C4.value,'OH- (M)':C5.value, 'I (M)':C6.value,'OCl- (M)':C7,'NH4+ (M)':C8,'H2CO3 (M)':C9,'HCO3 (M)':C10,'CO3 (M)':C11,'Alk (M)':C12}    

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)'])

标签： pythongekko

Answer 1

我想我找到了问题所在。我将 IMODE=4 更改为 IMODE=7 并且它起作用了：