Can Python SciPy odeint be run in a while loop?

Question

I am wondering if I can use SciPy's odeint in a while loop.

The example I am using is from http://apmonitor.com/pdc/index.php/Main/PhysicsBasedModels. I adjusted the code so that it is only solving the energy balance. The adjusted code is below and it works:

Code 1

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# define mixing model
def vessel(x,t,q,qf,Tf):
    # Inputs (3):
    # qf  = Inlet Volumetric Flowrate (L/min)
    # q   = Outlet Volumetric Flowrate (L/min)
    # Tf  = Feed Temperature (K)

    # States (2):
    # Volume (L)
    V = x[0]
    # Temperature (K)
    T = x[1]

    # Mass balance: volume derivative
    dVdt = qf - q

    # Energy balance: temperature derivative
    # Chain rule: d(V*T)/dt = T * dV/dt + V * dT/dt
    dTdt = (qf*Tf - q*T)/V - (T*dVdt/V)

    # Return derivatives
    return [dVdt,dTdt]

# Initial Conditions for the States
V0 = 1.0
T0 = 350.0
y0 = [V0,T0]

# Time Interval (min)
t = np.linspace(0,10,100)

# Inlet Volumetric Flowrate (L/min)
qf = np.ones(len(t))* 5.2
qf[50:] = 5.1

# Outlet Volumetric Flowrate (L/min)
q = np.ones(len(t))*5.0

# Feed Concentration (mol/L)
Caf = np.ones(len(t))*1.0
Caf[30:] = 0.5

# Feed Temperature (K)
Tf = np.ones(len(t))*300.0
Tf[70:] = 325.0

# Storage for results
V  = np.ones(len(t))*V0
T  = np.ones(len(t))*T0

# Loop through each time step
for i in range(len(t)-1):
    # Simulate
    inputs = (q[i],qf[i],Tf[i])
    ts = [t[i],t[i+1]]
    y = odeint(vessel,y0,ts,args=inputs)
    # Store results
    V[i+1]  = y[-1][0]
    T[i+1]  = y[-1][1]
    # Adjust initial condition for next loop
    y0 = y[-1]

# Plot the inputs and results
plt.figure()

plt.subplot(2,2,1)
plt.plot(t,qf,'b--',linewidth=3)
plt.plot(t,q,'b:',linewidth=3)
plt.ylabel('Flow Rates (L/min)')
plt.legend(['Inlet','Outlet'],loc='best')

plt.subplot(2,2,2)
plt.plot(t,Tf,'k--',linewidth=3)
plt.ylabel('Tf (K)')
plt.legend(['Feed Temperature'],loc='best')
plt.xlabel('Time (min)')

plt.subplot(2,2,3)
plt.plot(t,V,'b-',linewidth=3)
plt.ylabel('Volume (L)')
plt.legend(['Volume'],loc='best')

plt.subplot(2,2,4)
plt.plot(t,T,'k-',linewidth=3)
plt.ylabel('T (K)')
plt.legend(['Temperature'],loc='best')
plt.xlabel('Time (min)')

plt.show()

I would like the above code to work in a while loop and run every second using ts. My attempt is below:

Code 2

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from time import time, sleep, strftime, localtime

# define mixing model
def vessel(x,t,q,qf,Tf):
    # Inputs (3):
    # qf  = Inlet Volumetric Flowrate (L/min)
    # q   = Outlet Volumetric Flowrate (L/min)
    # Tf  = Feed Temperature (K)

    # States (2):
    # Volume (L)
    V = x[0]
    # Temperature (K)
    T = x[1]

    # Mass balance: volume derivative
    dVdt = qf - q

    # Energy balance: temperature derivative
    # Chain rule: d(V*T)/dt = T * dV/dt + V * dT/dt
    dTdt = (qf*Tf - q*T)/V - (T*dVdt/V)

    # Return derivatives
    return [dVdt,dTdt]

# Initial Conditions for the States
V0 = 1.0
T0 = 350.0
y0 = [V0,T0]

# Inlet Volumetric Flowrate (L/min)
qf = 5.2

# Outlet Volumetric Flowrate (L/min)
q = 5.0

# Feed Temperature (K)
Tf = 300.0

while 1:
    # Initial conditions
    V  = V0
    T  = T0
    
    t = strftime("%H:%M:%S", localtime()) #current local time
    ts = sum(int(x) * 60 ** i for i,x in enumerate(reversed(t.split(":")))) #current local time in seconds
    # Simulate
    inputs = (q,qf,Tf)
    ts = [t[i],t[i+1]]
    y = odeint(vessel,y0,ts,args=inputs)
    # Store results
    V[i+1]  = y[-1][0]
    T[i+1]  = y[-1][1]
    # Adjust initial condition for next loop
    y0 = y[-1]

My issue is in the while loop, I am confused on how to combine the while loop with the matrices used to obtain the values of the ODE. Every example I have come across so far using SciPy's odeint uses a set time interval like the for loop in Code 1. I am wondering if I need to use another tool like the GEKKO Optimization Suite (https://gekko.readthedocs.io/en/latest/), or maybe solve my equations a different way if I want to use a while loop.

Let me know if more detail is needed, or anything needs clarification.

Answer 1

Based on the additional comments, this modified for loop should solve the problem.

def get_t_local():
    t = strftime("%H:%M:%S", localtime()) #current local time
    return sum(int(x) * 60 ** i for i,x in enumerate(reversed(t.split(":")))) #current local time in seconds

times = [get_t_local()] ## starting time point
V, T = [V0], [T0]

inputs = (q,qf,Tf)
while 1:
    sleep(1)
    t_new = get_t_local() # get new timepoint

    # solve differential equation, take final result only
    y0 = odeint(vessel, y0, [times[-1], t_new], args = inputs) [-1]
    
    # Store results
    times.append(t_new)
    V.append(y0[0])
    T.append(y0[1])
    
    print(t_new, y0)