I am estimating some parameters with mode 5 in GEKKO and then passing these parameters to a simulation model (mode 4), but I don't get the same results for my variables. Is this caused by accumulated error due to the fact that the values of the parameters might get truncated or is there another issue that could be causing this? Should I get the same values for my variables? The variables are bounded (in case that affects the result).
Here is the code:
from gekko import GEKKO
def run_model_fixed(days, population, case, k_val, b_val, u0_val,
sigma_val,
Kmax0, a_val, c_val):
list_x =[]
list_u =[]
list_Kmax =[]
list_b =[]
error=[]
der=[]
list_s=[]
for i in range(len(days)):
list_xi=[]
list_ui=[]
list_bi=[]
list_Ki=[]
list_si=[]
list_der=[]
for j in range(len(days[i])):
m = GEKKO(remote=False)
m.time= days[i][j]
x_data= population[i][j]
x = m.CV(value=x_data, lb=0); x.FSTATUS = 1 # fit to measurement
x.SPLO = 0
sigma= m.Param(sigma_val)
s= m.Param(c_val)
d = m.Param(c_val)
k = m.FV(value=0.1, lb= 0, ub=100); k.STATUS=1
b = m.Param(b_val)
r = m.FV(value=0.5, lb= a_val, ub=1); r.STATUS=1
step = [0 if z<0 else 1 for z in m.time]
m_param = m.Param(step)
u = m.Var(value=u0_val, lb=0, ub=1)
m.free(u)
c= m.Param(c_val)
Kmax= m.Param(Kmax0)
if case == 'case4':
m.Equations([x.dt()== x*(r*(1-u**2)*(1-x/(Kmax))-m_param/(k+b*u)-d), u.dt() == sigma*(-2*u*r*(1-x/(Kmax))+(b*m_param)/(b*u+k)**2)])
m.options.IMODE = 5 # dynamic estimation
m.options.NODES = 5 # collocation nodes
m.options.EV_TYPE = 1 # linear error (2 for squared)
m.options.MAX_ITER=15000
m.solve(disp=False, debug=False) # display solver output
list_xi.append(x.value)
list_ui.append(u.value)
list_der.append(Kmax.value[0])
list_Ki.append(r.value[0])
list_bi.append(k.value[0])
list_si.append(b.value[0])
list_x.append(list_xi)
list_Kmax.append(list_Ki)
list_u.append(list_ui)
list_b.append(list_bi)
list_s.append(list_si)
der.append(list_der)
return list_x, list_u, list_Kmax, error, list_b, list_s, der
def run_model_sim(days, population, case, k_val, b_val, u0_val, sigma_val, Kmax0, a_val, c_val):
list_x =[]
list_u =[]
list_Kmax =[]
list_b =[]
error=[]
der=[]
list_s=[]
for i in range(len(days)):
list_xi=[]
list_ui=[]
list_bi=[]
list_Ki=[]
list_si=[]
list_der=[]
for j in range(len(days[i])):
#try:
m = GEKKO(remote=False)
m.time= days[i][j]
x_data= population[i][j]
x = m.Var(value=x_data[0], lb=0)
sigma= m.Param(sigma_val)
d = m.Param(c_val)
k = m.Param(k_val)
b = m.Param(b_val)
r = m.Param(a_val)
step = [0 if z<0 else 1 for z in m.time]
m_param = m.Param(step)
u = m.Var(value=u0_val, lb=0, ub=1)
#m.free(u)
a = m.Param(a_val)
c= m.Param(c_val)
Kmax= m.Param(Kmax0)
dx = m.Var(value = 0.0)
m.Equation(dx==x.dt())
if case == 'case4':
m.Equations([x.dt()== x*(r*(1-u**2)*(1-x/(Kmax))-m_param/(k+b*u)-d), u.dt() == sigma*(-2*u*r*(1-x/(Kmax))+(b*m_param)/(b*u+k)**2)])
m.options.IMODE = 5
#m.options.MAX_ITER=15000
#m.options.SOLVER=1
m.solve(disp=False, GUI=False)
list_xi.append(x.value)
list_ui.append(u.value)
#list_der.append(Kmax.value[0])
list_Ki.append(m_param.value)
#list_bi.append(k.value[0])
#list_si.append(sigma.value[0])
list_x.append(list_xi)
list_Kmax.append(list_Ki)
list_u.append(list_ui)
list_b.append(list_bi)
list_s.append(list_si)
der.append(list_der)
return list_x, list_u, list_Kmax, error, list_b, list_s, m.options.OBJFCNVAL
k0,b0,group, case0, u0, sigma0, K0, a0, c0 = (100, 100, 'Size1, Inc', 'case4', 0.1, 0.01, 0.1, 0, 0.01)
scaled_days[i][j]=[-7, 43, 104]
scaled_pop[i][j]=[0.00048029248653781915, 0.0005998737292033572, 0.0005998737292033572]
list_x, list_u, list_Kmax, error, list_b, list_s, der =run_model_fixed(days=[[scaled_days[i][j]]], population= [[scaled_pop[i][j]]],case=case0, k_val=k0, b_val=b0, u0_val=u0, sigma_val=sigma0, Kmax0=K0, a_val=a0, c_val=c0)
list_x2, list_u2, list_Kmax2, error2, list_b2, list_s2, der2 =run_model_sim(days=[[scaled_days[i][j]]], population= [[scaled_pop[i][j]]],case=case0, k_val=list_b[0][0], b_val=b0, u0_val=list_u[0][0], sigma_val=sigma0, Kmax0=K0, a_val=list_Kmax[0][0], c_val=c0)
If I check list_x and list_x2 are not the same