Pack voltage of 1s4p battery starkly different from cell voltages

[MustafaDragonborne]

I tried to simulate a battery pack using the Thevenin model with 4 cells in parallel in a 1s4p configuration and also managed to add in some slight soc imbalance between the cells. A drive cycle was used for the simulation. I notice however that the pack terminal voltage is very different from the cell terminal voltages with large errors between the two even though the two must be more or less the same since all cells are in parallel. The pack voltage also goes very frequently outside the allowed voltage limits of the cells (black dotted lines). Any suggestions here are highly appreciated.

EDIT: The last plot in the voltage subplots should be titled as 'Difference between pack voltage and mean cell voltage'

The code for the simulation is as follows:

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

import sys
sys.path.append('/home/mustafa/Documents/ECM_parameter_sets')
import liTec_Prototyp_40Ah_Zelle1

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib tk
from sklearn.metrics import mean_absolute_error, mean_squared_error

import pybamm
from pybamm import ParameterValues
pybamm.set_logging_level('NOTICE')
import liionpack as lp

curr = pd.read_csv('test_current_profile.csv')
curr['timestamp'] = pd.to_datetime(curr['timestamp'], unit='s')

# Create current values for the ECM model with LiTec cell having 40Ah capacity
curr_ecm = curr.copy()
curr_ecm['current'] = curr_ecm['current'] / 7.5  # 300/40 = 7.5 scaling factor
f'Max current value = {max(curr_ecm["current"])}'
f'Max C-rate = {max(curr_ecm["current"])/40}'  # cell nominal capacity = 40Ah
drive_cycle_ecm = curr_ecm.to_numpy(dtype='float64')[:,(0,2)]
Np = 4  # number of cells in parallel
Ns = 1  # number of cells in series
drive_cycle_ecm[:,1] = -drive_cycle_ecm[:,1] * Np

# set the parameter values according to chemistry
parameter_values_litec = liTec_Prototyp_40Ah_Zelle1.get_parameter_values()
# remove unwanted keys
unwanted_keys = 'chemistry'
del parameter_values_litec[unwanted_keys]
parameter_values_litec = ParameterValues(parameter_values_litec)
parameter_values_litec.update({'Initial SoC': '[input]'})
inputs = {'Initial SoC': np.linspace(0.5, 0.6, Np*Ns)}

# create the experiment
experiment = pybamm.Experiment(
    [("Run dc (A)")] * 1,
    drive_cycles={"dc": drive_cycle_ecm},
    period='1 second'
)
def mysim_func(parameter_values):
    # create the model
    model = pybamm.equivalent_circuit.Thevenin(options={'number of rc elements':2})

    # Set up solver and simulation
    solver = pybamm.CasadiSolver(mode="fast")
    sim = pybamm.Simulation(
        model=model,
        # experiment=experiment,
        parameter_values=parameter_values,
        solver=solver
    )
    return sim

netlist = lp.setup_circuit(
    Np=Np, Ns=Ns, Ri=1e-2, Rc=1e-2, Rb=1e-4, Rt=1e-5, I=0
)

output = lp.solve(
    netlist=netlist,
    sim_func=mysim_func,
    parameter_values=parameter_values_litec,
    experiment=experiment,
    inputs=inputs,
)