Estimation parameters using GEKKO (a second-order differential equation)

2
votes

I am using a second-order parallel model:

enter image description here

It has four parameters (initial values Cf0 and Cs0, kf and ks) to be estimated using GEKKO and six experiments. However I am not achieving a good fit, only experiments 1 and 4 are moving near the experiment data. Maybe I have some error in the code. Despite this, any guide-lines or link to get more information on controlling boundaries and the estimation method would be very helpful.

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

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

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

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

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

1 Answers

0
votes

The STATUS=1 needs to be before the m.solve() to make kf and ks adjustable by the solver.

if True:
    kf.STATUS=1
    ks.STATUS=1

The additional inequality constraint aren't needed because they already have lower bounds set elsewhere.

m.Equation(ks>0)
m.Equation(kf>0)
m.Equation(cnf01>0)
m.Equation(cns01>0)

The initial conditions are calculated by including fixed_initial=False.

cnf1=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns1=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf2=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns2=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf3=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns3=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf4=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns4=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf5=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns5=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cnf6=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns6=m.Var(value=1.3,lb=0.0,ub=10,fixed_initial=False)

One way to shorten the model is to use Arrays.

cnf1,cnf2,cnf3,cnf4,cnf5,cnf6 = m.Array(m.Var,6,\
            value=1.3,lb=0.0,ub=10,fixed_initial=False
cns1,cns2,cns3,cns4,cns5,cns6 = m.Array(m.Var,6,\
            value=1.3,lb=0.0,ub=10,fixed_initial=False

The IPOPT solver fails to converge after 250 iterations. Switching to the APOPT solver gives a successful solution.

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

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

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

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

tm6 = [0,0.08333,0.5,1,4.6667]
mca6 = [18.90, 15.8, 11.70,15.5,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,lb=0.01,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,cnf2,cnf3,cnf4,cnf5,cnf6 = m.Array(m.Var,6,\
            value=1.3,lb=0.0,ub=10,fixed_initial=False)
cns1,cns2,cns3,cns4,cns5,cns6 = m.Array(m.Var,6,\
            value=1.3,lb=0.0,ub=10,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,'b',label='Predicted c4')
plt.plot(m.time,c5.value,'r',label='Predicted c5')
plt.plot(m.time,c6.value,'g',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,'bo',label='Meas c4')
plt.plot(tm5,mca5,'ro',label='Meas c5')
plt.plot(tm6,mca6,'go',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