pythonscipyodeint可以在while循环中运行吗？

Question

我想知道是否可以在while循环中使用SciPy的odeint

我正在使用的示例来自http://apmonitor.com/pdc/index.php/Main/PhysicsBasedModels。我调整了代码，使其只解决了能量平衡问题。调整后的代码如下所示，且有效：

代码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()

我希望上面的代码在while循环中工作，并使用ts每秒运行一次。我的尝试如下：

代码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]

我的问题是在while循环中，我对如何将while循环与用于获取ODE值的矩阵组合感到困惑。到目前为止，我在使用SciPy的odeint时遇到的每个示例都使用一个设置的时间间隔，如代码1中的for循环。我想知道我是否需要使用另一种工具，比如GEKKO优化套件（https://gekko.readthedocs.io/en/latest/），或者如果我想使用while循环，也许可以用另一种方法来解决我的方程

如果需要更多的细节，或者有什么需要澄清的，请告诉我

Answer 1

根据附加注释，这个修改的for loop应该可以解决这个问题

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)