Linear Interpolation in the ODE (GEKKO)

Question

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:

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.

Answer 1

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).