Export Pybamm equations to Matlab / Simulink via Casadi: Temperature input not working

I would like to export the Pybamm equations as Casadi objects and solve them in Matlab. The code attached is mostly based on this example GitHub - FaradayInstitution/pybamm_simulink_example: An example of how to run pybamm in simulink. While my code works for current as the single input (marked as Version 1 in the Code), when adding the temperature input (Version2 in code) I get an error. Could somebody please help me out?

In the original example the temperature input is handled a little differently with an external model. However, it seems that this method is no longer available for more recent versions of Pybamm (24.9.0 used here).

To test the example:

  1. first run the python file Casadi_ModelCreation_example_SET_T.py to export the casadi objects

  2. run test_run_pybamm_casadi_SET_T.m which will call the function run_pybamm_casadi_example_SET_T.m for a few steps.

NOTE: Please specify the path / directory for your casadi exports and the casadi installation in run_pybamm_casadi_example_SET_T.m line 7 and 9 to make the code work!

# Casadi_ModelCreation_example_SET_T.py

import pybamm
import casadi as ca
import numpy as np
import time
import pickle
from scipy import io
import os
import shutil


options = {"thermal": "isothermal"} # lumped
model = pybamm.lithium_ion.DFN(options=options)

parameter_values = pybamm.ParameterValues("Chen2020") #('ORegan2022')#
T = 273.15 # in Kelvin = 0°C
parameter_values.update({"Ambient temperature [K]": T}) 
parameter_values.update({"Initial temperature [K]": T})
parameter_values["Lower voltage cut-off [V]"] = 2.0

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Version_1: Does work  
parameter_values.update({
        "Current function [A]": "[input]",
        #"Ambient temperature [K]": "[input]",   
    })

# Version_2: Does not work  
parameter_values.update({
        "Current function [A]": "[input]",
        "Ambient temperature [K]": "[input]",   
    })

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I_app = -5

solver = pybamm.CasadiSolver(mode="fast")

SoC_begin = 0.1 # 0%


T_in=273.15

sim = pybamm.Simulation(model, parameter_values=parameter_values, solver=solver)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Version_1: Does work  
inputs = {"Current function [A]": I_app,
           #"Ambient temperature [K]": T_in
           }

# Version_2: Does not work  
inputs = {"Current function [A]": I_app,
           "Ambient temperature [K]": T_in}

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


t_eval = np.linspace(0, 1, 2)

sim.solve(t_eval=t_eval, inputs=inputs,initial_soc=SoC_begin)
s=sim.solution
y0 = sim.solution.y.full()[:, 0] 

current_file_dir = os.path.dirname(os.path.abspath(__file__))
temp_dir = os.path.join(current_file_dir, 'casadi_export_SET_T')

shutil.rmtree(temp_dir, ignore_errors=True)
os.mkdir(temp_dir)
io.savemat(os.path.join(temp_dir, 'y0_SET_T.mat'), {'y0':y0})

casadi_integrator = solver.create_integrator(sim.built_model, inputs=inputs, t_eval=t_eval)
ci_path = os.path.join(temp_dir, 'integrator_SET_T.casadi')
casadi_integrator.save(ci_path)
variable_names = ['Terminal voltage [V]',
                  'Negative electrode surface potential difference at separator interface [V]',
                  'Volume-averaged cell temperature [C]',
                  'Current [A]'
                    ]


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Version_1: Does work  
ipo = ["Current function [A]"] # input parameter order

# Version_2: Does not work 
ipo = ["Current function [A]", "Ambient temperature [K]"] # input parameter order
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

casadi_objs = sim.built_model.export_casadi_objects(variable_names=variable_names,input_parameter_order=ipo)
variables = casadi_objs['variables']
t, x, z, p = casadi_objs["t"], casadi_objs["x"], casadi_objs["z"], casadi_objs["inputs"]
#rhs = casadi_objs["rhs"]
variables_stacked = ca.vertcat(*variables.values())
variables_fn = ca.Function("variables", [t, x, z, p], [variables_stacked])
v_path = os.path.join(temp_dir, 'variables_SET_T.casadi')
variables_fn.save(v_path)
print('Casadi objects re-generated \n', os.listdir(temp_dir))


% test_run_pybamm_casadi_SET_T.m 

clc

currents = [-5, -6, -7, -5, -6]; % current ampere

Ts = [25, 26, 27, 28, 26]+273.15;  % T /K 
total_steps = 5;

t_start = tic;

results_in =[];
obj_in=[];

for i = 1:total_steps
    % Perform a single time step simulation
    [results_out,obj_out]= run_pybamm_casadi_example_SET_T(results_in,obj_in,currents(i), Ts(i));
    obj_in = obj_out; 
    results_in=results_out;
end

t_end = toc(t_start);

print('done')
%run_pybamm_casadi_example_SET_T.m

function [results,obj] = run_pybamm_casadi_example_SET_T(results,obj,I_in,T_in)

% add casadi fodler to path 
addpath('C:\ CHANGE THIS \Casadi') % path of casasi installation

casadi_export_dir = 'C:\ CHANGE THIS  \casadi_export_SET_T'; % path of casadi objects
import casadi.*

if isempty(results) 

    % initialization
    tmp = load(fullfile(casadi_export_dir, 'y0_SET_T.mat'));
    obj.f = Function.load(fullfile(casadi_export_dir, 'integrator_SET_T.casadi'));
    obj.y0_1 = tmp.y0(1:862);
    obj.y0_2 = tmp.y0(862+1:end);

    obj.variables = Function.load(fullfile(casadi_export_dir, 'variables_SET_T.casadi'));
    
    
    
    %% Version1: only current input: Does work (when inputs created accordingly in python file)
    inputs = [I_in T_in]; 
    
    % Version2: only current input: Does  NOT work (when inputs created accordingly in python file)
    inputs = [I_in]; 

    % solve Model
    [yt_1,yt_2] = obj.f(obj.y0_1,obj.y0_2,horzcat(inputs,1e-6),0,0,0,0);   %t_min=1e-6

    % set new initial state for next step as last solved state
    obj.y0_1= yt_1(:,end);
    obj.y0_2= yt_2(:,end);
    % calc outputs
    outputs = double(full(obj.variables(0, obj.y0_1, obj.y0_2, inputs(1:end-1))));

    % initilaize and empty struct for the results
    results=struct();
    results.v = outputs(1);
    results.phi_an = outputs(2);
    results.T=outputs(3);
    results.I=outputs(4);
else    
    
    %% solve Model
    inputs = [I_in,1e-6]';

    [yt_1,yt_2] = obj.f(obj.y0_1,obj.y0_2, inputs,0,0,0,0);

    % set new initial state for next step as last solved state
    
    obj.y0_1= yt_1(:,end);
    obj.y0_2= yt_2(:,end);

    % calc outputs
    outputs = double(full(obj.variables(0, obj.y0_1, obj.y0_2, inputs(1:end-1))));
 
    results.v(end+1) = outputs(1);
    results.phi_an(end+1) = outputs(2);
    results.T(end+1)=outputs(3);
    results.I(end+1)=outputs(4);

end   


I’m interested in the same topic and have the same issue, could you share your experience with me if you have solved the problem or have been contacted by the PyBaMM team?