How to implement rectangular pulses (discontinuities) on ODE right-hand side?

In addition to this question, I would like to have several jumps during simulation time. As @LutzL proposed, I tried to perform one simulation per phase in a loop, and use the Connection method to stitch them together, with the final states of Phase 1 == initial states of Phase 2 etc. But I get an error, saying

Exception: @error: Model Expression *** Error in syntax of function string: Missing operator

Position: 6
11.55,11.55,11.55,11.55,11.55 ?

The feeding times and rates are to be understood as: [t_start, t_end], so there is two feeding events (starting at t = 2 and t = 3.1, respectively). one takes 0.7 and the other 0.2 (days). So there is five phases (numOfPhases): Before first feeding, first feeding, between first and second feeding, second feeding, after second feeding.

Here is my code:

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from gekko import GEKKO
# Construct input schedule
feedingTimes = np.array([2, 2.7, 3.1, 3.3])
feedingRates = np.array([60, 0, 430, 0])
tf      = 5
# Number of phases
numOfPhases = len(feedingTimes) + 1
# Construct time vectors (unequally spaced for now)
timeVector = np.zeros((numOfPhases, 10))
timeVector[0,:] = np.linspace(0, feedingTimes[0], 10)
for tI in np.arange(1,numOfPhases-1):
    timeVector[tI,:] = np.linspace(feedingTimes[tI-1], feedingTimes[tI], 10)
timeVector[-1,:] = np.linspace(feedingTimes[-1], tf, 10) 

#Initialize Model
m = GEKKO(remote=False)

q_in    = [m.Var(0.0, lb=-2, ub=500, fixed_initial=False) for i in range(numOfPhases)]

# (Example of) same parameter for each phase
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)
# Variables (one version of x for each phase)
X       = [m.Var(value = 11.55) for i in range(numOfPhases)]
Y       = [m.Var(value = 11.55*0.2) for i in range(numOfPhases)]

rho_1   = [m.Intermediate(k_1*X) for i in range(numOfPhases)]
q_prod  = [m.Intermediate(0.52*f_1*X) for i in range(numOfPhases)]

# Equations (different for each phase)
for pI in range(numOfPhases):
    m.time = timeVector[pI,:]
    m.Equations([X[pI].dt() == q_in[pI]/V_liq*(X_in - X[pI]) - rho_1[pI], \
                 Y[pI].dt() == q_in[pI]/V_liq*(Y_in - Y[pI])])

# Connect phases together at endpoints
for pI in range(numOfPhases-1):
    m.Connection(X[pI+1], X[pI], 1, len(m.time)-1, 1, 3)
    m.Connection(Y[pI+1], Y[pI], 1, len(m.time)-1, 1, 3)

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

plt.plot(m.time, q_in.value, label='q_in')
plt.plot(m.time, X.value, label='X')
plt.xlabel('time')
plt.ylabel('X')
plt.grid()
plt.legend(loc='best')
plt.show()

Solution

The Intermediates rho_1 and q_prod are not defined correctly because X is a variable list, not a single variable.

rho_1   = [m.Intermediate(k_1*X[i]) for i in range(numOfPhases)]
q_prod  = [m.Intermediate(0.52*f_1*X[i]) for i in range(numOfPhases)]

However, I don't think you need an array of X values for your problem because you want the final condition of one segment to equal the initial condition of the next segment and your equations do not change from one segment to the next. Multiple definitions of m.time are not allowed. You only define time once for each solve and each segment must have the same time. There is a way to time scale each segment so that it is variable but this is a more advanced topic. Also, the timeVector must not have two sequential values that are equal (must be increasing). I added 1e-5 to the start of each feeding phase.

timeVector = np.zeros((numOfPhases, 10))
timeVector[0,:] = np.linspace(0, feedingTimes[0], 10)
for tI in np.arange(1,numOfPhases-1):
    timeVector[tI,:] = np.linspace(feedingTimes[tI-1]+1e-5, feedingTimes[tI], 10)
timeVector[-1,:] = np.linspace(feedingTimes[-1]+1e-5, tf, 10)

The timeVector is a (5,10) numpy array. You need to flatten it into a 1D array for Gekko. I made this modification with the numpy.reshape function.

m.time = np.reshape(timeVector, -1)

If you need different equations for different time segments then you may need to define an array of X values with different equations. Because you use the same equations for each segment, I recommend a sequential approach instead as shown below.

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from gekko import GEKKO
# Construct input schedule
feedingTimes = np.array([2, 2.7, 3.1, 3.3])
feedingRates = np.array([60, 0, 430, 0])
tf      = 5
# Number of phases
numOfPhases = len(feedingTimes) + 1
# Construct time vectors (unequally spaced for now)
timeVector = np.zeros((numOfPhases, 10))
timeVector[0,:] = np.linspace(0, feedingTimes[0], 10)
for tI in np.arange(1,numOfPhases-1):
    timeVector[tI,:] = np.linspace(feedingTimes[tI-1]+1e-5, feedingTimes[tI], 10)
timeVector[-1,:] = np.linspace(feedingTimes[-1]+1e-5, tf, 10)

feedVector = np.zeros(tf*10)
for i in range(len(feedingRates)):
    feedVector[(i+1)*10:(i+2)*10] = feedingRates[i]

#Initialize Model
m = GEKKO(remote=False)
m.time = np.reshape(timeVector, -1)

q_in    = m.Param(value=feedVector)

# (Example of) same parameter for each phase
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)
# Variables (one version of x for each phase)
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)

# Equations (different for each phase)
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=True)

plt.plot(m.time, q_in.value, label='q_in')
plt.plot(m.time, X.value, label='X')
plt.xlabel('time')
plt.grid()
plt.legend(loc='best')
plt.show()

If you need to eventually switch to optimization / control and want only one value of feedrate for each segment then I recommend the parameter MV_STEP_HOR as described in the documentation.