Do I need to use Gekko m.connection in simulation mode to solve my code?

user_14965584 0
Viewed 69

To solve the question How to apply m.connection from Gekko using arrays?

I tried to use a new strategy:

I started testing if I can do a simulation before the estimation because I found some values to use as parameters. To simplify I simulate one data only and I assume ammonia is the only reactive reagent affecting the decay of the variable C2 (NH2Cl) by assuming parameters kdoc1 and kdoc2, doc1(t=0) and doc2(t=0) equal to zero, so r15 and r16 are also zero. The results here make sense because the decay is only due to ammonia. To get a higher decay and a better approximation to the data (squares in red), r15 and r16 should not be zero, indicating the contributions of organic matter in the decay of C2 (NH2Cl), too.

Here is the graph result:

enter image description here

Then my idea was to do the simulation adding values to the parameters (the values were validated). Even with DOF=0 I got an error 1. Is it possible I need to use m.connection to inform doc10 is doc1 at t=0 even in simulation mode? Or need to identify the variables as SV. or FV.? I try to apply but it did not work. Thank you in advance.

Here is my code assuming values for the parameters kdoc1, kdoc2, doc10 (doc1(t=0)) and 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,pos1=1,pos2=1,node1=1,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,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)'])
gekko python
1 Answers

To answer the question on how to set the initial conditions for dynamic, the value is the correct way to specify the initial condition in simulation mode such as x=m.Var(value=0).

doc1 = m.Var(value=doc10)

If the initial condition is calculated then set fixed_initial=False with IMODE=5 to allow degrees of freedom.

m.options.IMODE=5
doc1 = m.Var(value=doc10,fixed_initial=False)

The doc10 value gives an initial guess that is then changed by the optimizer to minimize the objective function such as the squared error fit to data.

Using IMODE=7 produces an equivalent result to IMODE=4 but solves the dynamic equations sequentially instead of simultaneously. It is often easier to solve.

user_23669410
Related