Different results giving same inputs

[MarcoZangrandi]

Hello everyone, I am trying to simulate a 3-cell parallel of the Chen2020 type. I wanted to pay attention to the first step.
I built a function that calculates the distribution of the initial currents in the parallel and gives the values you see in the first code, i.e. [[3.2500733514010998, 0.21684771789384336, -1.466921069294921]] respectively for cell 1, cell 2 and cell 3.
Then I simulated the first step for the three cells that have different initial SOCs.

With regard to the bulk open-circuit voltages, the results are as follows:

• Step 0: V_start = array([[4.2 , 4.09666099, 4.0420787 ]]) # [Cell1, Cell2, Cell3]
• Step 1: V_last = array([[4.18843279, 4.09553281, 4.03736792]]) # [Cell1, Cell2, Cell3]

Since it didn't make sense to me that cell 3, with a negative current discharged (voltage dropping) I tried simulating the same thing myself in another code (the second one you can find).
Actually providing the same current the voltage for the first step of cell 3 comes out like this:
V_bulk = array([4.0420787, 4.0424485]) # [Cell 3]
and correctly charges.

I don't understand why by providing the same values the simulations give two different results.
Can anyone answer me? Thank you

First code

import numpy as np
import math
import matplotlib.pyplot as plt
import pybamm
import warnings
import itertools
import time
from gen_initial_currents_2 import initial_current_guessing

# Cells configuration
np = 3  # cells in parallel
ns = 1  # cells in series
N_batteries = np.arange(0, np * ns) 
N_batteries = N_batteries.reshape(ns, np)
t_final = 3 # time in [h]
t_final = t_final*3600 # time in [s]
SOC_limit_low = 0.009 # lower limit for SOC
SOC_limit_high = 1 # upper limit for SOC

# Inputs
Ibatt = 2  # parallel group input current in [A]
SOC_start = [ [1, 0.9, 0.8] ] # initial cells SOC

dt = 5  # step time in [s] 

# Model definition
options = {"thermal": "lumped"}  # thermal model definition
model = pybamm.lithium_ion.SPMe(options=options) # electrochemical model definition
geometry = pybamm.battery_geometry(form_factor="cylindrical")  # geometry definition
param_sets = [["Chen2020", "Chen2020", "Chen2020"]] # cell sets definition 
param = [[0 for _ in range(np)] for _ in range(ns)]
sim = [[0 for _ in range(np)] for _ in range(ns)]

# Simulation objects creation (in this case it will create np simulation objects)
for cell_groups, param_groups in enumerate(param_sets):
    for cell_elem, param_elem in enumerate(param_groups):
        param[cell_groups][cell_elem] = pybamm.ParameterValues(param_elem)
        param[cell_groups][cell_elem] = pybamm.ParameterValues.set_initial_stoichiometries(param[cell_groups[cell_elem], initial_value=SOC_start[cell_groups][cell_elem])
        param[cell_groups][cell_elem]["Current function [A]"] = "[input]"
        sim[cell_groups][cell_elem] = pybamm.Simulation(model, geometry=geometry, parameter_values=param[cell_groups][cell_elem])
        cell_elem = 0
    cell_groups = 0

# Initial current guessing
initial_currents = initial_current_guessing(Ibatt, np, ns, dt, SOC_start, options, model, geometry, param, sim)
# the function "initial_currents_guessing" provides these values [[3.2500733514010998, 0.21684771789384336, -1.466921069294921]]

n = len(N_batteries)
m = len(N_batteries[0])
V_start = np.zeros_like(N_batteries, dtype=float)
Cap_start = np.zeros_like(N_batteries, dtype=float)
dSOC = np.zeros_like(N_batteries, dtype=float)
V_limit_high = np.zeros_like(N_batteries, dtype=float)
V_limit_low = np.zeros_like(N_batteries, dtype=float)
V_term_step0 = np.zeros_like(N_batteries, dtype=float)
V_term_step1 = np.zeros_like(N_batteries, dtype=float)
V_last = np.zeros_like(N_batteries, dtype=float)
SOC_last = np.zeros_like(N_batteries, dtype=float)
Ri_step0 = np.zeros_like(N_batteries, dtype=float)
Ri_step1 = np.zeros_like(N_batteries, dtype=float)

# First step
for i in range(n):
    for j in range(m):
        sim[i][j].step(dt=dt, npts=2, inputs={"Current function [A]": initial_currents[i][j]})
        V_start[i][j] = sim[i][j].solution["Bulk open-circuit voltage [V]"].entries[0]
        Cap_start[i][j] = param[i][j]['Nominal cell capacity [A.h]']
        dsoc = (initial_currents[i][j]*dt)/(3600*param[i][j]['Nominal cell capacity [A.h]'])
        dSOC[i][j] = dsoc # array that contains the SOC variations
        V_limit_high[i][j] = param[i][j]["Upper voltage cut-off [V]"]
        V_limit_low[i][j] = param[i][j]["Lower voltage cut-off [V]"]
        V_last[i][j] = sim[i][j].solution["Bulk open-circuit voltage [V]"].entries[-1]
        V_term_step0[i][j] = sim[i][j].solution["Terminal voltage [V]"].entries[0]
        V_term_step1[i][j] = sim[i][j].solution["Terminal voltage [V]"].entries[-1]
        Ri_step0[i][j] = abs(sim[i][j].solution["Local ECM resistance [Ohm]"].entries[0])
        Ri_step1[i][j] = abs(sim[i][j].solution["Local ECM resistance [Ohm]"].entries[-1])
        j = 0
    i = 0

Second code

from scipy.interpolate import interp1d
import numpy as np
import math
import matplotlib.pyplot as plt
import pybamm
import warnings
import itertools
import time

# Model definition
options = {"thermal": "lumped"}  
model = pybamm.lithium_ion.SPMe(options=options) 
geometry = pybamm.battery_geometry(form_factor="cylindrical") 
param = pybamm.ParameterValues("Chen2020")
pybamm.ParameterValues.set_initial_stoichiometries(param, initial_value=0.8)
param["Current function [A]"] = "[input]"
sim = pybamm.Simulation(model, geometry=geometry, parameter_values=param)
sim.step(dt=5, npts=2, inputs={"Current function [A]": -1.466921069294921})
V_bulk = sim.solution["Bulk open-circuit voltage [V]"].entries

[TomTranter]

The biggest problem was this line

param[cell_groups][cell_elem] = pybamm.ParameterValues.set_initial_stoichiometries(param[cell_groups[cell_elem], initial_value=SOC_start[cell_groups][cell_elem])

which should be

param[cell_groups][cell_elem].set_initial_stoichiometries(SOC_start[cell_groups][cell_elem])

Another problem was using np for number of batteries in parallel which replaces the import of numpy. It took me a while to figure out what you are doing with number of batteries. I would suggest wrapping the solving bits into a function that accepts a simulation and a current then returns the simulation or solution.

This seems to work now

import numpy as np
import pybamm


# Cells configuration
Np = 3  # cells in parallel
Ns = 1  # cells in series
N_batteries = int(Np * Ns)
t_final = 3 # time in [h]
t_final = t_final*3600 # time in [s]
SOC_limit_low = 0.009 # lower limit for SOC
SOC_limit_high = 1 # upper limit for SOC

# Inputs
Ibatt = 2  # parallel group input current in [A]
SOC_start = np.array([ [1, 0.9, 0.8] ]) # initial cells SOC

dt = 5  # step time in [s] 

# Model definition
options = {"thermal": "lumped"}  # thermal model definition
model = pybamm.lithium_ion.SPMe(options=options) # electrochemical model definition
geometry = pybamm.battery_geometry(form_factor="cylindrical")  # geometry definition
param_sets = [["Chen2020", "Chen2020", "Chen2020"]] # cell sets definition 
param = [[0 for _ in range(Np)] for _ in range(Ns)]
sim = [[0 for _ in range(Np)] for _ in range(Ns)]

# Simulation objects creation (in this case it will create np simulation objects)
for cell_groups, param_groups in enumerate(param_sets):
    for cell_elem, param_elem in enumerate(param_groups):
        param[cell_groups][cell_elem] = pybamm.ParameterValues(param_elem)
        param[cell_groups][cell_elem].set_initial_stoichiometries(SOC_start[cell_groups][cell_elem])
        param[cell_groups][cell_elem]["Current function [A]"] = "[input]"
        sim[cell_groups][cell_elem] = pybamm.Simulation(model, geometry=geometry, parameter_values=param[cell_groups][cell_elem])

# Initial current guessing
initial_currents = [[3.2500733514010998, 0.21684771789384336, -1.466921069294921]]


V_start = np.zeros_like(SOC_start, dtype=float)
Cap_start = np.zeros_like(SOC_start, dtype=float)
dSOC = np.zeros_like(SOC_start, dtype=float)
V_limit_high = np.zeros_like(SOC_start, dtype=float)
V_limit_low = np.zeros_like(SOC_start, dtype=float)
V_term_step0 = np.zeros_like(SOC_start, dtype=float)
V_term_step1 = np.zeros_like(SOC_start, dtype=float)
V_last = np.zeros_like(SOC_start, dtype=float)
SOC_last = np.zeros_like(SOC_start, dtype=float)
Ri_step0 = np.zeros_like(SOC_start, dtype=float)
Ri_step1 = np.zeros_like(SOC_start, dtype=float)

solutions = []
# First step
for i in range(Ns):
    for j in range(Np):
        sim[i][j].step(dt=dt, npts=2, inputs={"Current function [A]": initial_currents[i][j]})
        V_start[i][j] = sim[i][j].solution["Bulk open-circuit voltage [V]"].entries[0]
        Cap_start[i][j] = param[i][j]['Nominal cell capacity [A.h]']
        dsoc = (initial_currents[i][j]*dt)/(3600*param[i][j]['Nominal cell capacity [A.h]'])
        dSOC[i][j] = dsoc # array that contains the SOC variations
        V_limit_high[i][j] = param[i][j]["Upper voltage cut-off [V]"]
        V_limit_low[i][j] = param[i][j]["Lower voltage cut-off [V]"]
        V_last[i][j] = sim[i][j].solution["Bulk open-circuit voltage [V]"].entries[-1]
        V_term_step0[i][j] = sim[i][j].solution["Terminal voltage [V]"].entries[0]
        V_term_step1[i][j] = sim[i][j].solution["Terminal voltage [V]"].entries[-1]
        Ri_step0[i][j] = abs(sim[i][j].solution["Local ECM resistance [Ohm]"].entries[0])
        Ri_step1[i][j] = abs(sim[i][j].solution["Local ECM resistance [Ohm]"].entries[-1])
        solutions.append(sim[i][j].solution)

plot = pybamm.QuickPlot(solutions)
plot.dynamic_plot()

[TomTranter]

also you might want to check out liionpack for solving pack problems https://liionpack.readthedocs.io/en/latest/

[MarcoZangrandi]

Thank you very much for your help!!
I am working on a thesis that aims precisely at comparing the algorithm I created with liionpack, so I am consulting it.
Sorry about the variable np, I didn't think it would interfere with numpy, I had modified it quickly just to write the post as in my code many variables were in Italian.
Now it works correctly if I enter the initial currents as an array of numbers.
If I use my function "initial_current_guessing" it gives me problems. I have pasted below the function I used, following Valentin Sulzer's advice in #4078. I initially create graphs where on the x-axis I have the initial currents (step 0) and on the y-axis the initial resistances (step 0). Then I iterate using Kirchhoff's laws and applying Thevenin until convergence.
Could it be that by carrying out this function, having called the simulation object as in the main programme, this could give problems once I return to the main programme and so I end up with different values for that reason?

(Update: I have tried changing the name of the variables, including those of the simulation object within my function but it still gives the initial problem even though the values of the initial currents provided by my 'initial_current_guessing' function are the same as those used before.)

from scipy.interpolate import interp1d
import numpy as np
import math
import matplotlib.pyplot as plt
import pybamm
import warnings

def initial_current_guessing(Ibatt, Np, Ns, dt, SOC_start, options, model, geometry, param, sim):
 # Creation of I-R graphs for subsequent interpolation
    current_points = np.linspace(0.0001, 20, 100) 
    initial_resistances = np.zeros([len(current_points), Ns * Np])
    initial_vocv = np.zeros([Ns, Np])
    for i in range(Ns):
        for j in range(Np):
            for index_cp, cp in enumerate(current_points):
                sim[i][j].solve(t_eval=[0, dt], initial_soc=SOC_start[i][j], inputs={"Current function [A]": cp})
                k = i*Np + j
                initial_resistances[index_cp][k] = (sim[i][j].solution["Local ECM resistance [Ohm]"].entries[0])
    if i==len(range(Ns))-1:
        for i in range(Ns):
            for j in range(Np):
                initial_vocv[i][j] = sim[i][j].solution["Bulk open-circuit voltage [V]"].entries[0]    
                
    # Calculation of the current in the cells at the initial step
    initial_current_step = Ibatt/Np
    initial_currents = [[initial_current_step] * Np for _ in range(Ns)]
    iter_res = [[0] * Np for _ in range(Ns)]

    tol = np.ones(Ns)
    steps_count = 0

    for parallel_groups in range(Ns):
        
        initial_step_currents = initial_currents[parallel_groups]
        
        while tol[parallel_groups]>1e-15:
            iter_i = parallel_groups*Np
            # for the current values I have in the vector initial_step_currents calculates the associated resistance by interpolating
            for index_cell, value_i in enumerate(initial_step_currents): 
                value_i = abs(value_i)
                resistance_points = [element[iter_i] for element in initial_resistances] # takes into account one column at a time of the 'initial_resistances' matrix
                # Plots I-R graphs
                plt.plot(current_points, resistance_points)
                plt.scatter(value_i, iter_res[parallel_groups][index_cell], color='red')
                interp_R = interp1d(current_points, resistance_points, kind='cubic') # interpolating function creation
                iter_res[parallel_groups][index_cell] = interp_R(value_i) # acquisition of the initial resistance as the current changes
                print(iter_res[parallel_groups][index_cell])
                iter_i+=1
                
            counter_cells = 0
            currents_prev = list(initial_step_currents) 
            while counter_cells < num_celle_p:
                # Thevenin
                r1_init = iter_res[parallel_groups][counter_cells]  # [Ohm]
                r_eq_init_inv = 0  # [Ohm]
                for counter_r_cells in range(len(iter_res[parallel_groups])):
                    if counter_r_cells != counter_cells:
                        r_eq_init_inv = r_eq_init_inv + 1 / iter_res[parallel_groups][counter_r_cells]
                r_eq_init = 1/r_eq_init_inv
                Vocv1 = initial_vocv[parallel_groups][counter_cells]
                i2 = 0.0  # initialization
                for counter_r_cells in range(len(iter_res[parallel_groups])):
                    if counter_r_cells != counter_cells:
                        i2 += initial_vocv[parallel_groups][counter_r_cells] / iter_res[parallel_groups][counter_r_cells]
                Vocv2 = i2/r_eq_init_inv
                
                # Kirchhoff's Law
                initial_currents[parallel_groups][counter_cells] = (Vocv1 - Vocv2 + Ibatt*r_eq_init)/(r1_init+r_eq_init)
                counter_cells += 1
            # Check on the tolerance 
            tol[parallel_groups] = min(abs(np.array(initial_currents[parallel_groups])-np.array(currents_prev))) # difference between current and previous iteration
            steps_count+=1
            
    return initial_currents

[TomTranter]

I'm afraid I've helped all I can