Parameters estimation regarding constraints for initial values are not working with m.connection or fixed_initial, any suggestion on that?

user_14965584 0
Viewed 80

Following the suggestion from Rate coefficients and initial concentration estimation using gekko - advice about model structure I adapted the code for two data sets and it ran (I am using COLDSTART=2).

enter image description here

However, I still have these doubts and it would be great to have some suggestions:

  1. I want to estimate kdoc1 and kdoc2 and the initial values (t=0) of DOC1 and DOC2 which are DOC10 and DOC20, regarding DOC10[1] = DOC10[0]*0.5 and DOC20[1] = DOC20[0]*0.5. I am expecting to get the relation for t=0 and after that DOC10 and DOC20 will decay, but it is not happening, values are keeping constant.

  2. I am not using m.connection because I got this error:

p13 not found in results file
p14 not found in results file
p15 not found in results file
p16 not found in results file

  1. and neither using fixed_initial =False because it gives infeasibilities

  2. Another issue is that I need to comment the variables initialization to achieve a result.

  3. The estimated values are equal to the initial values proposed. It seems the program is not fitting and not giving an estimation value.

Here is the code:

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

#data set 1
t_data1 = [0,0.08333,0.5,1,4,22.61667]
x_data1 = [0,4.91e-5,4.57e-5,4.74e-5,4.17e-5,2.76e-5]
x_data1mgl = [0,3.48,3.24,3.36,2.96,1.96]

#data set 2
t_data2 = [0,0.08333,0.5,1,4,22.8167]
x_data2 = [0,5.92e-5,5.7e-5,5.64e-5,5.30e-5,4.6e-5]
x_data2mgl = [0,4.2,4.04,4,3.76,2.88]
      
#combine data and time in the same dataframe
data1 = pd.DataFrame({'time':t_data1,'x1':x_data1})
data2 = pd.DataFrame({'time':t_data2,'x2':x_data2})
data1.set_index('time', inplace=True)
data2.set_index('time', inplace=True)


# merge dataframes
data = data1.join(data4, how='outer')

# simulation time points
dftp = pd.DataFrame({'time':np.linspace(0,25,51)})
dftp.set_index('time', inplace=True)

# merge dataframes
df = data.join(dftp,how='outer')

# get True (1) or False (0) for measurement
#df['meas'] = (df['x'].values==df['x'].values).astype(int)
z1 = (df['x1']==df['x1']).astype(int)
z2 = (df['x2']==df['x2']).astype(int)


# replace NaN with zeros
df0 = df.fillna(value=0)

#Estimator Model
m = GEKKO()#remote=True)
#m = GEKKO(remote=False) # remote=False to produce local folder with results

m.time = df.index.values

# measurements - Gekko arrays to store measured data
xm = m.Array(m.Param,2)
xm[0].value = df0['x1'].values
xm[1].value = df0['x2'].values

hocl_init_val=m.Array(m.Param,2)
hocl_init_val[0].value= 8.01e-5
hocl_init_val[1].value= 8.01e-5

nh3_init_val=m.Array(m.Param,2)
nh3_init_val[0].value= 1.37e-3 
nh3_init_val[1].value= 6.82e-4

h_init_val=m.Array(m.Param,2)
h_init_val[0].value= 6.31e-8 
h_init_val[1].value= 2e-8

#h2co3_init_val=m.Array(m.Param,2)
#h2co3_init_val[0].value= 5.39e-4
#h2co3_init_val[1].value= 1.88e-4

hco3_init_val=m.Array(m.Param,2)
hco3_init_val[0].value= 3.46e-3
hco3_init_val[1].value= 3.80e-3

#co32_init_val=m.Array(m.Param,2)
#co32_init_val[0].value= 2.16e-6
#co32_init_val[1].value= 7.5e-6

alk_init_val=m.Array(m.Param,2)
alk_init_val[0].value= 3.46e-3
alk_init_val[1].value= 3.82e-3

#adjustable parameters

x1=3.1684e-5
x2=3.1308e-5

DOC10=m.Array(m.FV,2)
DOC10[0].value= x1
DOC10[1].value= x1*0.5

DOC20=m.Array(m.FV,2)
DOC20[0].value= x2
DOC20[1].value= x2*0.5


kdoc1 = m.FV(1e6,lb=0.01,ub=1e7) #m.Const(1.2334e6)
kdoc2 = m.FV(1e9,lb=0.01,ub=1e10) #m.Const(3.8809e9)
#DOC10 = m.FV(1e-5,lb=1e-7,ub=1e-2)
#DOC20 = m.FV(1e-5,lb=1e-7,ub=1e-2)


#variables initialization
# hocl=m.Array(m.Var,2,lb=0,ub=1e-4)
# nh3=m.Array(m.Var,2,lb=0,ub=1.5e-4)
# nh2cl=m.Array(m.Var,2,lb=0,ub=1e-4)
# C2m=m.Array(m.Var,2)
# cC2=m.Array(m.Var,2)
# nhcl2=m.Array(m.Var,2)
# h=m.Array(m.Param,2)
# #oh=m.Array(m.Var,2)
# I=m.Array(m.Var,2)
# ocl=m.Array(m.Var,2)
# nh4=m.Array(m.Var,2)
# #h2co3=m.Array(m.Var,2)
# hco3=m.Array(m.Var,2)
# #co32=m.Array(m.Var,2)
# alk=m.Array(m.Var,2)
# DOC1=m.Array(m.SV,2)
# DOC2=m.Array(m.SV,2)

#to support intermediate
co32=[None]*2
h2co3=[None]*2

oh=[None]*2

r1=[None]*2
r2=[None]*2
r3=[None]*2
r4=[None]*2
r5=[None]*2
r6=[None]*2
r7=[None]*2
r8=[None]*2
r9=[None]*2
r10=[None]*2
r11=[None]*2
r12=[None]*2
r13=[None]*2
r14=[None]*2
r15=[None]*2
r16=[None]*2

#Define GEKKO variables that determine if time 
#point contains data to be used in regression
zm = m.Array(m.Param,2)
zm[0].value=z1
zm[1].value=z2

# fit to measurement
x=m.Array(m.Param,2) 
x[0].value=x_data1[0]
x[1].value=x_data4[0]

k5=m.Array(m.Var,2)

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)

for i in range(2):

    hocl[i] = m.Var(value=hocl_init_val[i],lb=0,ub=1e-2)
    nh3[i] = m.Var(value=nh3_init_val[i],lb=0,ub=1.5e-3)
    nh2cl[i] = m.Var(value=x[i],lb=0,ub=3e-4)
    
    C2m[i] = m.Param(xm[i],lb=0)
    meas = m.Param(zm[i])
    m.Minimize(meas*((nh2cl[i]-C2m[i]))**2)
    
    nhcl2[i] = m.Var(value=0)
    
    h[i] = m.Param(value=h_init_val[i])
    
    I[i] = m.Var(value=0)
    ocl[i]= m.Var(value=0)
    nh4[i] = m.Var(value=0)
    #h2co3[i] = m.Var(value=h2co3_init_val[i])
    hco3[i] = m.Var(value=hco3_init_val[i])
    #co32[i] = m.Var(value=co32_init_val[i])
    alk[i] = m.Var(value=alk_init_val[i])
    DOC1[i] = m.SV(value=DOC10[i])#,fixed_initial=False)
    DOC2[i] = m.SV(value=DOC20[i])#,fixed_initial=False)
    cC2[i] = m.Var(value=0)
    
    oh[i] = m.Intermediate(1e-14/h_init_val[i])

    r1[i] = m.Intermediate(k1 * hocl[i] * nh3[i])
    r2[i] = m.Intermediate(k2 * nh2cl[i])
    r3[i] = m.Intermediate(k3 * hocl[i] * nh2cl[i])
    r4[i] = m.Intermediate(k4 * nh3[i])
    r5[i] = m.Intermediate(k5[i] * nh2cl[i] * nh2cl[i])
    r6[i] = m.Intermediate(k6 * nhcl2[i] * nh3[i]* h[i])
    r7[i] = m.Intermediate(k7 * nhcl2[i] * oh[i])
    r8[i] = m.Intermediate(k8 * I[i] * nhcl2[i])
    r9[i] = m.Intermediate(k9 * I[i] * nh2cl[i])
    r10[i] = m.Intermediate(k10 * nh2cl[i] * nhcl2[i])
    
    r11[i] = m.Intermediate(k11*hocl[i])
    r12[i] = m.Intermediate(k12*h[i]*ocl[i])
    r13[i] = m.Intermediate(k13*nh4[i])
    r14[i] = m.Intermediate(k14*h[i]*nh3[i])
    
    r15[i] = m.Intermediate(kdoc1*DOC1[i]*nh2cl[i])
    r16[i] = m.Intermediate(kdoc2*DOC2[i]*hocl[i])

    co32[i] = m.Intermediate(5.01e-11*hco3[i]/h[i])
    h2co3[i] = m.Intermediate(hco3[i]*h[i]/5.01e-7)
    
    t = m.Param(value=m.time)
    m.Equation(hocl[i].dt()== -r1[i] + r2[i] - r3[i] + r4[i] + r8[i] - r11[i] + r12[i] - r16[i])
    m.Equation(nh3[i].dt()== -r1[i] + r2[i] + r5[i] - r6[i] + r13[i] - r14[i])# + r11[i])
    m.Equation(nh2cl[i].dt()== r1[i] - r2[i] - r3[i] + r4[i] - r5[i] + r6[i] - r9[i] - r10[i] - r15[i])
    m.Equation(nhcl2[i].dt()== r3[i] - r4[i] + r5[i] - r6[i] - r7[i] - r8[i] - r10[i])
    #m.Equation(h[i].dt()== 0)
    #m.Equation(oh[i] == 1e-14/h[i])
    m.Equation(I[i].dt()== r7[i] - r8[i] - r9[i])
    #m.Equation(ocl[i] == (3.16e-8*hocl[i])/h[i])
    #m.Equation(nh4[i] == nh3[i]*h[i]/5.01e-10)
    m.Equation(ocl[i].dt()==r11[i]-r12[i])
    m.Equation(nh4[i].dt()==-r13[i]+r14[i])
    #m.Equation(co32[i] == (5.01e-11*hco3[i])/h[i])
    #m.Equation(h2co3[i] == (hco3[i])*h[i]/5.01e-7)
    m.Equation(hco3[i] == alk[i] - 2*5.01e-11*hco3[i]/h[i] - oh[i] + h[i])
    m.Equation(alk[i].dt()== 0)
    
    m.Equation(DOC1[i].dt()== -r15[i])
    m.Equation(DOC2[i].dt()== -r16[i])
    
    m.Equation(cC2[i] == 51500*nh2cl[i])
    
    m.Equation(k5[i] == 2.5e7*h[i]+4e4*(hco3[i]*h[i]/5.01e-7)+800*hco3[i])

    #m.Connection(DOC1[i],DOC10[i],pos1=1,pos2=1,node1=1,node2=1)
    #m.Connection(DOC2[i],DOC20[i],pos1=1,pos2=1,node1=1,node2=1)
    
#Application options
m.options.SOLVER = 3 #IPOPT solver
m.options.EV_TYPE = 2 #absolute error
#m.options.NODES = 3 #collocation nodes (2,5)
#m.options.RTOL = 1E-6

if True:
    kdoc1.STATUS=1
    kdoc2.STATUS=1
    DOC10[i].STATUS=1
    DOC20[i].STATUS=1

#m.options.IMODE = 7
#m.options.RTOL = 1E-6
#m.solve(dis=True)

m.options.IMODE = 5 #Dynamic Simultaneous - estimation = MHE
m.options.COLDSTART=2
m.options.TIME_SHIFT=0

m.open_folder()
m.solve(disp=True)

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

print('Solution')
print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
print('DOC10_w1 = ' + str(DOC10[0].value[0]))
print('DOC20_w1 = ' + str(DOC20[0].value[0]))
print('DOC10_w2 = ' + str(DOC10[1].value[1]))
print('DOC20_w2 = ' + str(DOC20[1].value[1]))

plt.figure(1,figsize=(8,5))
plt.subplot(2,1,1)
plt.plot(m.time,nh2cl[0],'b.--',label='Predicted 1')
plt.plot(m.time,nh2cl[1],'r.--',label='Predicted 2')

plt.plot(t_data1,x_data1,'bx',label='Meas 1')
plt.plot(t_data2,x_data2,'rx',label='Meas 2')
#plt.legend(loc='best')
plt.ylabel('mol/L')
plt.legend(loc='center right', bbox_to_anchor=(1.45, 0.5), ncol=2) #shadow=True,

plt.subplot(2,1,2)
plt.plot(m.time,cC2[0].value,'b',label ='C2_M1')
plt.plot(m.time,cC2[1].value,'r',label ='C2_M2')
plt.plot(t_data1,x_data1mgl,'bx',label='Meas 1')
plt.plot(t_data2,x_data2mgl,'rx',label='Meas 2')

plt.ylabel('mg Cl2/L')
plt.xlabel('time (h)')
plt.legend(loc='center right', bbox_to_anchor=(1.4, 0.5),  ncol=2) #shadow=True,


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

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

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

plt.subplot(4,3,4)
plt.plot(m.time,h[i].value,label ='h')
plt.legend()

plt.subplot(4,3,5)
plt.plot(m.time,oh[i].value,label ='oh')
plt.legend()

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

plt.subplot(4,3,7)
plt.plot(m.time,ocl[i].value,label ='ocl')
plt.legend()

plt.subplot(4,3,8)
plt.plot(m.time,nh4[i].value,label ='nh4')
plt.legend()

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

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

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

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

plt.show()
gekko
1 Answers

Here are two ideas:

  1. Try a range of fixed initial conditions DOC10 and DOC20 and then re-optimize with every combination. This can be parallelized with multithreading so that multiple optimization problems run simultaneously to speed up the solution.
import numpy as np
from gekko import GEKKO

# Optimize at mesh points
x = np.arange(20.0, 30.0, 2.0)
y = np.arange(30.0, 50.0, 4.0)
amg, bmg = np.meshgrid(x, y)

# Initialize results array
obj = np.empty_like(amg)

m = GEKKO(remote=False)
objective = float('NaN')

a,b = m.Array(m.FV,2)

# model variables, equations, objective
x1 = m.Var(1,lb=1,ub=5)
x2 = m.Var(5,lb=1,ub=5)
x3 = m.Var(5,lb=1,ub=5)
x4 = m.Var(1,lb=1,ub=5)
m.Equation(x1*x2*x3*x4>=a)
m.Equation(x1**2+x2**2+x3**2+x4**2==b)
m.Minimize(x1*x4*(x1+x2+x3)+x3)
m.options.SOLVER = 1 # APOPT solver

# Calculate obj at all meshgrid points
for i in range(amg.shape[0]):
    for j in range(bmg.shape[1]):
        a.MEAS = amg[i,j]
        b.MEAS = bmg[i,j]

        m.solve(disp=False)

        obj[i,j] = m.options.OBJFCNVAL
        print(amg[i,j],bmg[i,j],obj[i,j])

# plot 3D figure of results
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(amg, bmg, obj, \
       rstride=1, cstride=1, cmap=cm.coolwarm, \
       vmin = 10, vmax = 25, linewidth=0, antialiased=False)
ax.set_xlabel('a')
ax.set_ylabel('b')
ax.set_zlabel('obj')
plt.show()

range of optimization values

  1. If the optimizer needs to chose the values, another option is to write out the discretization equations with orthogonal collocation on finite elements. This gives much more control with IMODE=3 because all of the equations and variables are individually defined for each time point. Here is a Gekko model with IMODE=4:
from gekko import GEKKO
m = GEKKO(remote=False) # create GEKKO model
y = m.Var(5.0,name='y') # create GEKKO variable
m.Equation(y.dt()==-y)  # create GEKKO equation
m.time = [0,1,2,3]      # time points
m.options.IMODE = 4     # simulation mode
m.options.NODES = 4     # set NODES=4
m.options.CSV_WRITE=2   # write results_all.json
                        #   with internal nodes
m.solve(disp=False)     # solve

and the equivalent model with IMODE=3:

import numpy as np
N = np.array([[0.436,-0.281, 0.121], \
              [0.614, 0.064, 0.0461], \
              [0.603, 0.230, 0.167]])

time = np.array([0.0, \
                 0.5-np.sqrt(5)/10.0, \
                 0.5+np.sqrt(5)/10.0, \
                 1.0])

from gekko import GEKKO
m = GEKKO()
y0 = 5
y  = m.Array(m.Var,3,value=0)
dy = m.Array(m.Var,3,value=0)

Ndy = np.dot(N,dy)
m.Equations([Ndy[i]==(y[i]-y0) for i in range(3)])
m.Equations([dy[i]+y[i]==0 for i in range(3)])
m.solve(disp=False)
print(y)

This would cause a major rewrite of the model but may be necessary if more control is needed. Additional information on Orthogonal Collocation on Finite Elements is in the Dynamic Optimization course.

user_23669410
Related