Dynamic optimization with linear interpolation for some parameters in GEKKO
This is a follow-up question from my previous post on linear interpolation while solving ODEs in GEKKO. Basically, I have a parameter k as a linear function of time, k(t), like:
and I want to include this parameter to be automatically interpolated during ODE integration.
Below is the code for running the simulation:
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()
My follow-up question is that I now want to calculate k(t), i.e., optimization, to achieve some objective function.
So, I modified the above code to change k from Param to MV and switched to IMODE = 6.
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.MV(value=k,lb = 0, ub = 2,fixed_initial=False)
k.STATUS = 1
x = m.Var(value=0)
m.Equation(x.dt() == k)
m.options.IMODE = 6
m.Minimize(x)
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.show()
The thing is that the above code calculates k for 13 points corresponding to the length of m.time. However, I want the solver to calculate k for only 5 points corresponding to the original length of kv. The reason is that, in the actual problem, only 5 points are enough for optimziation, and so calculating 13 points is too expensive. This is a cartoon example, so the reason may not be obvious, but in the real problem, it is 15 points for kv vs 150 points for m.time (the ODE is very stiff), and so that is really important.
Another approach I tried is force k to have 5 points:
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 = m.MV(value=kv, lb = 0, ub = 2,fixed_initial=False)
k.STATUS = 1
x = m.Var(value=0)
m.Equation(x.dt() == k)
m.options.IMODE = 6
m.Minimize(x)
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.show()
But got an error: Data arrays must have the same length, and match time discretization in dynamic problems, which is understandable.
Really appreciate if someone could help solve this problem.
Thank you!
Set m.options.MV_STEP_HOR to adjust from the default value of 1 (make a move every step) to a higher value, such as 4 to make a manipulated variable move every 4 steps. To keep the same slope in that interval, define a new manipulated variable ks as the slope and k as a variable.
k = m.Var(lb=0, ub=2, fixed_initial=False)
ks = m.MV(value=0,lb=-10,ub=10); ks.STATUS=1
m.Equation(k.dt()==ks)
The lower lb and upper ub bounds on ks set the minimum and maximum rate of change of k. The objective function is modified to make the result more interesting. The blue dashed is the slope ks and the red line is k. The slope of k changes when ks changes every 4 time points. There are 10 time steps in this problem so the changes occur at t=0, t=0.4, t=0.8 with 0.1 time intervals.
m.Minimize((x-0.5)**2)
Here is the full script:
from gekko import GEKKO
import numpy as np
m = GEKKO()
m.time = np.linspace(0, 1, 11)
k = m.Var(lb=0, ub=2, fixed_initial=False)
ks = m.MV(value=0,lb=-10,ub=10); ks.STATUS=1
m.Equation(k.dt()==ks)
x = m.Var(value=0)
m.Equation(x.dt() == k)
m.options.IMODE = 6
m.options.MV_STEP_HOR=4
m.Minimize((x-0.5)**2)
m.solve(disp=False)
import matplotlib.pyplot as plt
plt.figure(figsize=(5,2.5))
plt.step(m.time,ks.value,'b:',label='ks')
plt.plot(m.time,k.value,'rx-',label='k')
plt.plot(m.time,x.value,'k-',label='x')
plt.plot(); plt.legend(); plt.grid()
plt.show()