Casadi integrator unable to perform a dynamic simulation in Pyomo

Question

I'm trying to perform a dynamic simulation of the filling process of a hydrogen storage tank in Pyomo. My code is as follows:

from pyomo.environ import *
from pyomo.dae import *
import numpy as np
from pem_simulation import T, p_H2_co, n_H2_co          # Outputs from another model which define the initial conditoons for some variables

# %% Model creation
sp = ConcreteModel()

sp.t_end = Param(initialize=(5000))
sp.t = ContinuousSet(bounds = (0, sp.t_end))

# %% Parameters

sp.V = Param(initialize=(0.411))                        # Vessel volume [m3] 
sp.R = Param(initialize=(8.314))                        # Ideal gas constant [J/(mol K)] 
sp.c_p = Param(initialize=(28.824))                     # Specific heat capacity at const. pressure [J/(mol K)] 
sp.c_v = Param(initialize=(sp.c_p - sp.R))              # Specific heat capacity at const. volume [J/(mol K)] 
sp.kappa = Param(initialize=sp.c_p/sp.c_v)              # Isentrope exponent [-]

# %% Variables

# Model variables

sp.T = Var(sp.t, initialize = 298.15)                             # Vessel temperature [K]
sp.p = Var(sp.t, initialize = 1e5)                                # Vessel pressure [Pa]
sp.n = Var(sp.t, initialize = sp.p[0]*sp.V / (sp.R*sp.T[0]))      # Amount of hydrogen in the vessel [mol]

sp.T_dot = DerivativeVar(sp.T)        
sp.n_dot = DerivativeVar(sp.n)

sp.P_ein = Var(sp.t, initialize = 0)                              # Compressor power consumption [W]
sp.p_el = Var(sp.t, initialize = p_H2_co[-1]*10**5)               # Elektorlyseur output pressure [Pa] 
sp.T_el = Var(sp.t, initialize = T[-1])                           # Elektrolyseur output temperature [K]
sp.n_e = Var(sp.t, initialize = n_H2_co[-1])                      # Hydrogen input amount [mol/s] 

# Auxiliary variables
sp.tau_e = Var(sp.t, initialize=(sp.n[0]/sp.n_e[0]))              # Inflow time parameter [s]
sp.T_0 = Var(sp.t, initialize=(sp.T[0]))                          # Initial temperature [K]
sp.T_t = Var(sp.t, initialize=(sp.R*sp.T[0] / sp.c_v))            # Temperature parameter for the energy balance [K]

# %% Constraint functions

def T_0f(sp, t):                                                                    # Function for T_0
    return sp.T_0[t] == T[0]                                                 
def T_tf(sp, t):                                                                    # Function for T_t
    return sp.T_t[t] == sp.R * sp.T[0] / sp.c_v

def iggf(sp, t):                                                                    # Ideal gas equation function
    return sp.p[t] * sp.V == sp.n[t] * sp.R * sp.T[t]

def n_inf(sp, t):                                                                   # Input hydrogen flow function
    return sp.n_e[t] == n_H2_co[-1]
def  p_inf(sp, t):                                                                  # Input pressure function
    return sp.p_el[t] == p_H2_co[-1]*10**5
def  T_inf(sp, t):                                                                  # Input temperature function
    return sp.T_el[t] == T[-1]

def tau_inf(sp, t):                                                                 # Function for tau_in
    return sp.tau_e[t] == sp.n[t]/sp.n_e[t]

def mb_inf(sp, t):                                                                  # Mass balance function
    return sp.n_dot[t] == sp.n_e[t] 
def eb_inf(sp, t):                                                                  # Energy balance function
    return sp.T_dot[t] == (sp.kappa * (sp.T_el[t] - sp.T_0[t]) - (sp.T[t] - sp.T_0[t]) + sp.T_t[t]) / sp.tau_e[t]

def P_inf(sp, t):                                                                   # Function for P_in
    return sp.P_ein[t] == sp.n_e[t] * sp.R * sp.T_el[t] * log(sp.p[t] / sp.p_el[t])

# %% Equations defined as constraints

sp.T_0c = Constraint(sp.t, rule = T_0f)
sp.T_tc = Constraint(sp.t, rule = T_tf)

sp.iggc = Constraint(sp.t, rule = iggf)

sp.tau_inc = Constraint(sp.t, rule = tau_inf)

sp.n_inc = Constraint(sp.t, rule = n_inf)
sp.p_inc = Constraint(sp.t, rule = p_inf)
sp.T_inc = Constraint(sp.t, rule = T_inf)

sp.mb_inc = Constraint(sp.t, rule = mb_inf)
sp.eb_inc = Constraint(sp.t, rule = eb_inf)

sp.P_inc = Constraint(sp.t, rule = P_inf)

# %% Simulation

sim = Simulator(sp, package = 'casadi')
tsim, profiles = sim.simulate(numpoints=100, integrator = 'collocation')

I need to post the whole code because I don't have an idea where exactly the problem might be. When I try to perform it, I receive the following warning:

WARNING("F_rootfinder:jac_f_z failed: NaN detected for output jac_o1_i1, at nonzero index 2 (row 5, col 0).") [D:\bld\casadi_1647512547064\work\casadi\core\oracle_function.cpp:265]

As a result of this, my solution matrix for the variable values at each time step contains mostly zeros plus the initial conditions for sp.T and sp.n repeated at each time step. I'm experiencing the same problem in another model as well.

I don't know what more to do in order to repair this, so I would appreciate any help. I'm using Python V. 3.10 and Anaconda Navigator (anaconda3).