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:
-
first run the python file Casadi_ModelCreation_example_SET_T.py to export the casadi objects
-
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