Incremental solving of GEKKO dynamic models

Question

I would like to use GEKKO to incrementally solve a dynamic model. I need to insert some parameters at different points in the simulation that are determined stochastically using an external model.

As a simple proof of concept, I tried adjusting the PID control example from the tutorial website to solve in 10 minute increments.

m = GEKKO()

tf = 40
timesteps = np.linspace(0,tf,2*tf+1)
step = np.zeros(2*tf+1)
step[3:40] = 2.0
step[40:]  = 5.0

# Controller model
Kc = 15.0                    # controller gain
tauI = 2.0                  # controller reset time
tauD = 1.0                  # derivative constant
OP_0 = m.Const(value=0.0)   # OP bias
OP = m.Var(value=0.0)       # controller output
PV = m.Var(value=0.0)       # process variable
SP = m.Param(value=step)    # set point
Intgl = m.Var(value=0.0)    # integral of the error
err = m.Intermediate(SP-PV) # set point error
m.Equation(Intgl.dt()==err) # integral of the error
m.Equation(OP == OP_0 + Kc*err + (Kc/tauI)*Intgl - PV.dt())

# Process model
Kp = 0.5                    # process gain
tauP = 10.0                 # process time constant
m.Equation(tauP*PV.dt() + PV == Kp*OP)

setpoint_final = []
pv_final = []
output_final = []
for i in range(8):
    SP.value = step[10*i:10*(i+1)]
    m.time = timesteps[10*i:10*(i+1)]
    m.options.IMODE=7
    m.solve(disp=False)
    setpoint_final.append(SP.value)
    pv_final.append(PV.value)
    output_final.append(OP.value)

plt.figure()
plt.subplot(2,1,1)
plt.plot(timesteps[:-1],np.concatenate(output_final),'b:',label='OP')
plt.ylabel('Output')
plt.legend()
plt.subplot(2,1,2)
plt.plot(timesteps[:-1],np.concatenate(setpoint_final),'k-',label='SP')
plt.plot(timesteps[:-1],np.concatenate(pv_final),'r--',label='PV')
plt.xlabel('Time (sec)')
plt.ylabel('Process')
plt.legend()
plt.show()

I should get something like this:

Instead it seems to reset at each interval:

Is there something I am missing or is this type of simulation not possible in GEKKO?

Answer 1

Try using m.options.TIME_SHIFT=10, m.options.IMODE=4, and change the way the values are stored.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO(remote=False)

tf = 40
timesteps = np.linspace(0,tf,2*tf+1)
step = np.zeros(2*tf+1)
step[3:40] = 2.0
step[40:]  = 5.0

# Controller model
Kc = 15.0                    # controller gain
tauI = 2.0                  # controller reset time
tauD = 1.0                  # derivative constant
OP_0 = m.Const(value=0.0)   # OP bias
OP = m.Var(value=0.0)       # controller output
PV = m.Var(value=0.0)       # process variable
SP = m.Param(value=step)    # set point
Intgl = m.Var(value=0.0)    # integral of the error
err = m.Intermediate(SP-PV) # set point error
m.Equation(Intgl.dt()==err) # integral of the error
m.Equation(OP == OP_0 + Kc*err + (Kc/tauI)*Intgl - PV.dt())

# Process model
Kp = 0.5                    # process gain
tauP = 10.0                 # process time constant
m.Equation(tauP*PV.dt() + PV == Kp*OP)

setpoint_final = []
pv_final = []
output_final = []
m.options.TIME_SHIFT=10
m.time = timesteps[0:10]
for i in range(8):
    SP.value = step[10*i:10*(i+1)]
    m.options.IMODE=4
    m.solve(disp=False)
    setpoint_final = np.concatenate((setpoint_final,SP.value))
    pv_final = np.concatenate((pv_final,PV.value))
    output_final = np.concatenate((output_final,OP.value))

plt.figure()
plt.subplot(2,1,1)
plt.plot(timesteps[:-1],output_final,'b:',label='OP')
plt.ylabel('Output')
plt.legend()
plt.subplot(2,1,2)
plt.plot(timesteps[:-1],setpoint_final,'k-',label='SP')
plt.plot(timesteps[:-1],pv_final,'r--',label='PV')
plt.xlabel('Time (sec)')
plt.ylabel('Process')
plt.legend()
plt.show()