How to define time-dependent, discrete parameter?

Question

Recently, I have built a small model with GEKKO. It contains a Parameter which actually changes with time. How can I implement that? I tried using if3, but it gives an error.

Here's the MWE:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Started on 10-08-2019

@author: winkmal
"""
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
#Initialize Model
m = GEKKO(remote=False)

# Parameters
k_1     = m.Param(value = 0.19)
f_1     = m.Param(value = 29.0)
V_liq   = m.Param(value = 159.0)

q_in    = m.Param(value = 2.5)
X_in    = m.Param(value = 271.77)
Y_in    = m.Param(value = 164.34)

X       = m.Var(value = 11.55)
Y       = m.Var(value = 11.55*0.2)
rho_1   = m.Intermediate(k_1*X)
q_prod  = m.Intermediate(0.52*f_1*X)

m.time  = np.arange(0,5,1/12)

m.Equations([X.dt() == q_in/V_liq*(X_in - X) - rho_1, \
             Y.dt() == q_in/V_liq*(Y_in - Y)]) 
#Dynamic simulation
m.options.IMODE = 4
m.solve(disp=False)

plt.plot(m.time, X.value)
plt.xlabel('time')
plt.ylabel('X')
plt.show()

I tried the following:

q_in    = m.if3(m.time - 2, 0, 2.5)

so that q_in would be 0 initially, and become 2.5 at time = 2. But I get the following error:

  File "/usr/local/lib/python3.7/site-packages/gekko/gekko.py", line 1838, in solve
    raise Exception(apm_error)

Exception:  @error: Equation Definition
 Equation without an equality (=) or inequality (>,<)
 (((1-int_v5))*([-2.-1.91666667-1.83333333-1.75-1.66666667-1.58333333
 STOPPING...

Do you have an idea how I can achieve this? Actually, this variable jumps several times between 0 and 60, and I have the time points available in a CSV file. Ideally, I could create a loop that would check at every iteration if it's time for q_in to change, and overwrite the current value accordingly.

Answer 1

You can read the input from a CSV and assign the time-varying values to q_in.value either during Parameter initialization (see Example #1) or else in a loop where the value changes each time integration interval (see Example #2). Examples 1 and 2 both produce the following result but Example 1 is faster.

Example 1 may also be faster with option m.options.IMODE=7 if you have a very long time horizon. IMODE=7 uses a sequential solution method instead of a simultaneous solution method.

Example 1

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
t = np.arange(0,5,1/12)
step = [0 if z<2 else 2.5 for z in t]
m = GEKKO(remote=False)
k_1     = m.Param(value = 0.19)
f_1     = m.Param(value = 29.0)
V_liq   = m.Param(value = 159.0)
q_in    = m.Param(value = step)
X_in    = m.Param(value = 271.77)
Y_in    = m.Param(value = 164.34)
X       = m.Var(value = 11.55)
Y       = m.Var(value = 11.55*0.2)
rho_1   = m.Intermediate(k_1*X)
q_prod  = m.Intermediate(0.52*f_1*X)
m.time  = t
m.Equations([X.dt() == q_in/V_liq*(X_in - X) - rho_1, \
             Y.dt() == q_in/V_liq*(Y_in - Y)]) 
m.options.IMODE = 4
m.solve(disp=False)
plt.plot(m.time,q_in.value,label=r'$q_{in}$')
plt.plot(m.time, X.value,label='X')
plt.plot(m.time, Y.value,label='Y')
plt.legend()
plt.xlabel('time')
plt.show()

Example 2

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
t = np.arange(0,5,1/12)
m = GEKKO(remote=False)
k_1     = m.Param(value = 0.19)
f_1     = m.Param(value = 29.0)
V_liq   = m.Param(value = 159.0)
q_in    = m.Param()
X_in    = m.Param(value = 271.77)
Y_in    = m.Param(value = 164.34)
X       = m.Var(value = 11.55)
Y       = m.Var(value = 11.55*0.2)
rho_1   = m.Intermediate(k_1*X)
q_prod  = m.Intermediate(0.52*f_1*X)
m.time  = [t[0],t[1]]
m.Equations([X.dt() == q_in/V_liq*(X_in - X) - rho_1, \
             Y.dt() == q_in/V_liq*(Y_in - Y)]) 
m.options.IMODE = 4
# store Xs and Ys for plotting
for i in range (1,len(t)):
    q_in.value = 0 if t[i]<2 else 2.5
    m.solve(disp=False)
    if i==1:
        Xs = [X.value[0]]
        Ys = [Y.value[0]]
    Xs.append(X.value[1])
    Ys.append(Y.value[1])
step = [0 if z<2 else 2.5 for z in t]
plt.plot(t,step,label=r'$q_{in}$')
plt.plot(t, Xs,label='X')
plt.plot(t, Ys,label='Y')
plt.legend()
plt.xlabel('time')
plt.show()

If you need to make q_in dependent on the value of some of your variables then you can use the m.if3 function. However, this is a more challenging problem to solve because the m.if3 function converts the problem into a Mixed Integer Nonlinear Programming form that may take longer to solve. Here is an example where q_in=0when X>8 and q_in=2.5 when X<=8. However, it didn't converge for me. I'm not sure why and I'd need to do some additional digging but I though you'd like to have it in case it works for you.

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=False)
k_1     = m.Param(value = 0.19)
f_1     = m.Param(value = 29.0)
V_liq   = m.Param(value = 159.0)
X_in    = m.Param(value = 271.77)
Y_in    = m.Param(value = 164.34)
X       = m.Var(value = 11.55,name='X')
Y       = m.Var(value = 11.55*0.2,name='Y')
rho_1   = m.Intermediate(k_1*X)
q_prod  = m.Intermediate(0.52*f_1*X)
q_in    = m.if3(8-X, 0.0, 2.5)
m.time  = np.arange(0,5,1/12)
m.Equations([X.dt() == q_in/V_liq*(X_in - X) - rho_1, \
             Y.dt() == q_in/V_liq*(Y_in - Y)]) 
m.options.IMODE = 6
m.options.SOLVER = 1
m.solve(disp=True)
plt.plot(m.time,q_in.value,label=r'$q_{in}$')
plt.plot(m.time, X.value,label='X')
plt.plot(m.time, Y.value,label='Y')
plt.legend()
plt.xlabel('time')
plt.show()

There are also a few other examples here on solving ODEs with time-varying inputs with Gekko.