Should I get the same results with the estimation mode and the simultaion mode in GEKKO?

Question

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

One other thing to check is that the solution parameters are the same between the simulator and optimizer. Values such as `m.options.NODES` can slightly adjust the solution values. You may also want to do a steady state initialization for both before starting. — John Hedengren, Dec 09 '20 at 21:14
Thanks, I added the same number of nodes and now the difference is much lower. With respect to the steady state initialization, I thought GEKKO automatically did that. How can I do this? — user606273, Dec 11 '20 at 08:57
For steady state initialization you can either (1) solve with `IMODE=1` first before solving in the other mode or (2) solve with `IMODE=4+` multiple times until the values reach steady state. — John Hedengren, Dec 11 '20 at 12:56
Thanks, and this is only to get the initial value of the variables? Meaning only x and u? And then if instead of simulating I am estimating parameters, do I also have to get the estimations with steady state and then use them to initialize values in the dynamic problem? Or only for simulation? — user606273, Dec 11 '20 at 14:13
Steady state initialization is typically only for simulation where you fix the `u` value and calculate the response `x`. — John Hedengren, Dec 11 '20 at 17:52

score 1 · Answer 1 · answered Dec 04 '20 at 16:00

As far as I understand, your results between estimation and simulation modes must be identical. You might want to check if you retrieved the estimated value correctly. The first value of your estimated variable or parameter array can remain as an initial guess value.

If you estimate a variable (m.Var) that can change throughout the simulation time, you can check if you used the whole array of your estimated variable when you evaluate the estimated parameter. However, if you estimate a parameter (m.Param), it stays at a single value throughout the whole simulation.

If you can post your code, I can give you a better answer.

I added a minimal reproducible example – user606273 Dec 04 '20 at 18:47 — user606273, Dec 04 '20 at 18:47

Should I get the same results with the estimation mode and the simultaion mode in GEKKO?

1 Answers1