Linear Interpolation in the ODE (GEKKO)

I have spent some time to solve this issue but still cannot find the solution. Basically, I have a parameter, say k, as a piecewise linear function of time, for example:

k time-varying

And so we can write k as an array vs time - k(t). The thing is I have an ODE system containing k, let say a simple one like: dx/dt = k(t), and I would like GEKKO to linear interpolate k from the above piecewise linear function during the ODE integration. This is basically a simulation problem.

After that, I will look at varying k(t) (i.e., optimal control using a piecewise linear control) to meet some objective function.

Is there any suggestion on how to do these things properly in GEKKO?

Thank you

I tried np.interp, but it did not accept GK variables. I also tried to write my own code for linear interpolation, but could not figure out how to extract the time point (m.time) during the integration.

Solution

Create the linear interpolation before creating the m.Param().

kt = np.array([0,0.2,0.45,0.75,1.0])
kv = np.array([0.5,1.0,0.2,1.2,1.0])
k = np.interp(m.time, kt, kv)
k = m.Param(value=k)

Here is an example solving at 0.1 time intervals from 0 to 1.

from gekko import GEKKO
import numpy as np

kt = np.array([0,0.2,0.45,0.75,1.0])
kv = np.array([0.5,1.0,0.2,1.2,1.0])

m = GEKKO()
m.time = np.linspace(0, 1, 11)

k = np.interp(m.time, kt, kv)
k = m.Param(value=k)

x = m.Var(value=0)
m.Equation(x.dt() == k)

m.options.IMODE = 4
m.solve(disp=False)

import matplotlib.pyplot as plt
plt.figure(figsize=(5,2.5))
plt.plot(kt,kv,'bo-',label='k data')
plt.plot(m.time,k.value,'rx-',label='k interp')
plt.plot(m.time,x.value,'k-',label='x')
plt.plot(); plt.legend(); plt.grid()
plt.tight_layout(); plt.savefig('p.png',dpi=300)
plt.show()

If there is interpolation error (as shown in the figure above), increase the number of time points or include the kt time points in m.time.

# include kt time points in m.time
m.time = np.concatenate((m.time, kt))
m.time = np.sort(m.time)
m.time = np.unique(m.time)

from gekko import GEKKO
import numpy as np

kt = np.array([0,0.2,0.45,0.75,1.0])
kv = np.array([0.5,1.0,0.2,1.2,1.0])

m = GEKKO()
m.time = np.linspace(0, 1, 11)
# include kt time points in m.time
m.time = np.concatenate((m.time, kt))
m.time = np.sort(m.time)
m.time = np.unique(m.time)

k = np.interp(m.time, kt, kv)
k = m.Param(value=k)

x = m.Var(value=0)
m.Equation(x.dt() == k)

m.options.IMODE = 4
m.solve(disp=False)

import matplotlib.pyplot as plt
plt.figure(figsize=(5,2.5))
plt.plot(kt,kv,'bo-',label='k data')
plt.plot(m.time,k.value,'rx-',label='k interp')
plt.plot(m.time,x.value,'k-',label='x')
plt.plot(); plt.legend(); plt.grid()
plt.tight_layout(); plt.savefig('p.png',dpi=300)
plt.show()

If k is calculated, set k=m.MV() with k.STATUS=1 for linear interpolation of the interior collocation time points. There is a global option m.options.MV_TYPE for either zero-order hold or first-order hold (linear interpolation).