Update current at each time step for a 2P configuration

[wigging]

I am trying to use the PyBaMM SPMe model in a parallel circuit. I created an example (see below) that represents two cells in parallel. To ensure the cells will behave differently, I changed the capacity of each cell using the Electrode height [m] parameter. The current in each branch is determined using a minimization function from SciPy. At each time step, I update the current applied to each SPMe model using the Current function [A] input. This updated branch current should result in the same voltage output from each cell. However, based on the plots shown below, it doesn't appear that the code is working correctly because the voltages are different. I have used this approach for solving branch currents with a parallel circuit of resistors without any issues. Maybe this approach does not work well when you represent the resistors as batteries? Or am I not using PyBaMM correctly to update the current applied to each cell? Any suggestions on how to run two SPMe models in a parallel (2P) configuration would be helpful.

import matplotlib.pyplot as plt
import numpy as np
import pybamm
from scipy.optimize import minimize

plt.style.use('default')

# Objective and constraint functions
# ----------------------------------------------------------------------------


def objfunc(x, rb):
    # Function to minimize voltage differences.
    v = x * rb
    v[1:] = v[1:] * -1
    return sum(v)**2


def con1(x, itot):
    # Constraint where sum of branch currents must equal total current.
    return sum(x) - itot


def con2(x, rb):
    # Constraint where all branch voltages must be equal.
    v = x * rb
    return np.diff(v)


# Parameters
# ----------------------------------------------------------------------------

dt = 20       # time step [s]
tf = 4000     # final time [s]
nb = 2        # number of parallel branches
itot = 1.0    # total current [A]

# Cell model
# ----------------------------------------------------------------------------

model = pybamm.lithium_ion.SPMe()

param1 = model.default_parameter_values
param1['Electrode height [m]'] = 0.134
param1["Current function [A]"] = "[input]"

param2 = model.default_parameter_values
param2['Electrode height [m]'] = 0.144
param2["Current function [A]"] = "[input]"

# Run simulation
# ----------------------------------------------------------------------------

# store cell 1 results
t1 = []
v1 = []
i1 = []

# store cell 2 results
t2 = []
v2 = []
i2 = []

# bounds for solving each branch current
bnds = [(0, itot) for _ in range(nb)]

# initial current in each branch [A]
ib0 = itot / nb
ib = [ib0, ib0]

# setup simulation objects for each cell
sim1 = pybamm.Simulation(model, parameter_values=param1)
sim2 = pybamm.Simulation(model, parameter_values=param2)

# simulation time [s]
time = 0

while time < tf:

    # solve a single time step
    sim1.step(dt, save=False, starting_solution=sim1.solution, inputs={'Current function [A]': ib[0]})
    sim2.step(dt, save=False, starting_solution=sim2.solution, inputs={'Current function [A]': ib[1]})

    # calculate branch resistances, R = V / I
    r1 = sim1.solution['Terminal voltage [V]'].entries[0] / sim1.solution['Current [A]'].entries[0]
    r2 = sim2.solution['Terminal voltage [V]'].entries[0] / sim2.solution['Current [A]'].entries[0]
    rb = [r1, r2]

    # constraints for current and voltage
    cons = [
        {'type': 'eq', 'fun': con1, 'args': (itot,)},
        {'type': 'eq', 'fun': con2, 'args': (rb,)}
    ]

    # determine corrected branch currents
    x0 = ib
    res = minimize(objfunc, x0, args=(rb,), method='SLSQP', bounds=bnds, constraints=cons)
    ib = res.x

    # update cell 1
    t1.append(sim1.solution['Time [s]'].entries[0])
    v1.append(sim1.solution['Terminal voltage [V]'].entries[0])
    i1.append(sim1.solution['Current [A]'].entries[0])

    # update cell 2
    t2.append(sim2.solution['Time [s]'].entries[0])
    v2.append(sim2.solution['Terminal voltage [V]'].entries[0])
    i2.append(sim2.solution['Current [A]'].entries[0])

    time += dt

# Plot
# ----------------------------------------------------------------------------

_, ax = plt.subplots(tight_layout=True)
ax.plot(t1, v1, label='cell 1')
ax.plot(t2, v2, label='cell 2')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Terminal voltage [V]')
ax.legend(loc='best')

_, ax = plt.subplots(tight_layout=True)
ax.plot(t1, i1, label='cell 1')
ax.plot(t2, i2, label='cell 2')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Current [A]')
ax.legend(loc='best')

plt.show()

[valentinsulzer]

Changing it to

# calculate branch resistances, R = (OCV - V) / I
# in pybamm this is the "Local ECM resistance [Ohm]" variable
r1 = sim1.solution["Local ECM resistance [Ohm]"].entries[0]
r2 = sim2.solution["Local ECM resistance [Ohm]"].entries[0]

makes it work but I can't get my head around why. @TomTranter has done some pack modeling and might be able to help.

Also your optimization setup seems weird. Having all the voltages be equal is both a constraint in con2 and the objective? And why not just define the last current in terms of all the others in order to satisfy con1 without needing an explicit constraint?

[TomTranter]

Yeah the ECM resistance is what you should use. It's like the internal resistance and needs to go in series with the voltage source. I'll be releasing my pack code next week

[wigging]

@tinosulzer Woo hoo! Your suggestion works. Thanks again for the help. I was not aware of the Local ECM resistance [Ohm] parameter. I'm still learning about all the features of PyBaMM. Also, I agree that the optimization procedure to solve for the branch currents could be improved. I just put something together that seemed to work for now.

@TomTranter Thank you for the comment. I meant to ask you earlier today about how you were solving for the currents but we ran out of time. We can have that discussion some other time. Until then, I put together this example to get more familiar with PyBaMM and to test out my approach.