Different Initial condition for the cells in pack

[ksnvikrant]

I tried to populate different initial electrolyte concentration across the cells in the pack. I have an array of initial electrolyte concentration (initial_eleConc) similar to total heat transfer coefficient (htc). The changes I made to utils.py, simulations.py, solver_utils.py are given below.

Changed line 220, and added line 229 in utils.py.

# -*- coding: utf-8 -*-
"""
Created on Thu Sep 23 10:33:13 2021

@author: Tom
"""
from scipy.interpolate import interp1d, interp2d
import numpy as np
import pandas as pd
import os
import liionpack
from skspatial.objects import Plane
from skspatial.objects import Points

ROOT_DIR = os.path.dirname(os.path.abspath(liionpack.__file__))
CIRCUIT_DIR = os.path.join(ROOT_DIR, "circuits")
DATA_DIR = os.path.join(ROOT_DIR, "data")
INIT_FUNCS = os.path.join(ROOT_DIR, "init_funcs")

def interp_current(df):
    r"""
    Returns an interpolation function for current w.r.t time

    Parameters
    ----------
    df : pandas.DataFrame or Dict
        Contains data for 'Time' and 'Cells Total Current' from which to
        construct an interpolant function

    Returns
    -------
    f : function
        interpolant function of total cell current with time.

    """
    t = df["Time"]
    I = df["Cells Total Current"]
    f = interp1d(t, I)
    return f

def _z_from_plane(X, Y, plane):
    r"""
    Given X and Y and a plane provide Z
    X - temperature
    Y - flow rate
    Z - heat transfer coefficient

    Parameters
    ----------
    X : float (array)
        x-coordinate.
    Y : float (array)
        z-coordinate.
    plane : skspatial.object.Plane
        plane returned from read_cfd_data.

    Returns
    -------
    z : float (array)
        z-coordinate.

    """
    a, b, c = plane.normal
    d = plane.point.dot(plane.normal)
    z = (d - a * X - b * Y) / c
    return z

def _fit_plane(xv, yv, dbatt):
    r"""
    Private method to fit plane to CFD data

    Parameters
    ----------
    xv : ndarray
        temperature meshgrid points.
    yv : ndarray
        flow_rate meshgrid points.
    dbatt : ndarray
        cfd data for heat transfer coefficient.

    Returns
    -------
    plane : skspatial.object.Plane
        htc varies linearly with temperature and flow rate so relationship
        describes a plane

    """
    nx, ny = xv.shape
    pts = []
    for i in range(nx):
        for j in range(ny):
            pts.append([xv[i, j], yv[i, j], dbatt[i, j]])

    points = Points(pts)
    plane = Plane.best_fit(points, tol=1e-6)
    return plane

def read_cfd_data(data_dir=None, filename="cfd_data.xlsx", fit="linear"):
    r"""
    A very bespoke function to read heat transfer coefficients from an excel
    file

    Parameters
    ----------
    data_dir : str, optional
        Path to data file. The default is None. If unspecified the module
        liionpack.DATA_DIR folder will be used
    filename : str, optional
        The default is 'cfd_data.xlsx'.
    fit : str
        options are 'linear' (default) and 'interpolated'.
    Returns
    -------
    funcs : list
        an interpolant is returned for each cell in the excel file.

    """
    if data_dir is None:
        data_dir = liionpack.DATA_DIR
    fpath = os.path.join(data_dir, filename)
    ncells = 32
    flow_bps = np.array(pd.read_excel(fpath, sheet_name="massflow_bps", header=None))
    temp_bps = np.array(pd.read_excel(fpath, sheet_name="temperature_bps", header=None))
    xv, yv = np.meshgrid(temp_bps, flow_bps)
    data = np.zeros([len(temp_bps), len(flow_bps), ncells])
    fits = []
    for i in range(ncells):
        data[:, :, i] = np.array(
            pd.read_excel(fpath, sheet_name="cell" + str(i + 1), header=None)
        )
        # funcs.append(interp2d(xv, yv, data[:, :, i], kind='linear'))
        if fit == "linear":
            fits.append(_fit_plane(xv, yv, data[:, :, i]))
        elif fit == "interpolated":
            fits.append(interp2d(xv, yv, data[:, :, i], kind="linear"))

    return data, xv, yv, fits

def get_linear_htc(planes, T, Q):
    r"""
    A very bespoke function that is called in the solve process to update the
    heat transfer coefficients for every battery - assuming linear relation
    between temperature, flow rate and heat transfer coefficient.

    Parameters
    ----------
    planes : list
        each element of the list is a plane equation describing linear relation
        between temperature, flow rate and heat transfer coefficient.
    T : float array
        The temperature of each battery.
    Q : float
        The flow rate for the system.

    Returns
    -------
    htc : float
        Heat transfer coefficient for each battery.

    """
    ncell = len(T)
    htc = np.zeros(ncell)
    for i in range(ncell):
        htc[i] = _z_from_plane(T[i], Q, planes[i])
    return htc

def get_interpolated_htc(funcs, T, Q):
    r"""
    A very bespoke function that is called in the solve process to update the
    heat transfer coefficients for every battery

    Parameters
    ----------
    funcs : list
        each element of the list is an interpolant function.
    T : float array
        The temperature of each battery.
    Q : float
        The flow rate for the system.

    Returns
    -------
    htc : float
        Heat transfer coefficient for each battery.

    """
    ncell = len(T)
    htc = np.zeros(ncell)
    for i in range(ncell):
        htc[i] = funcs[i](T[i], Q)
    return htc

def build_inputs_dict(I_batt, htc,initial_eleConc):
    r"""
    Function to convert inputs and external_variable arrays to list of dicts
    As expected by the casadi solver. These are then converted back for mapped
    solving but stored individually on each returned solution.
    Can probably remove this process later

    Parameters
    ----------
    I_batt : float array
        The input current for each battery.
    htc : float array
        the heat transfer coefficient for each battery.

    Returns
    -------
    inputs_dict : list
        each element of the list is an inputs dictionary corresponding to each
        battery.

    """
    inputs_dict = []
    for i in range(len(I_batt)):
        inputs_dict.append(
            {
                # 'Volume-averaged cell temperature': T_batt[i],
                "Current": I_batt[i],
                "Total heat transfer coefficient [W.m-2.K-1]": htc[i],
                "Initial concentration in electrolyte [mol.m-3]": initial_eleConc[i],
            }
        )
    return inputs_dict

Added line 73 in simulations.py.

# -*- coding: utf-8 -*-
"""
Created on Wed Sep 22 15:37:51 2021

@author: tom
"""

import pybamm

def _current_function(t):
    r"""
    Internal function to make current an input parameter

    Parameters
    ----------
    t : float
        Dummy time parameter.

    Returns
    -------
    TYPE
        DESCRIPTION.

    """
    return pybamm.InputParameter("Current")

def create_simulation(parameter_values=None, experiment=None, make_inputs=False):
    r"""
    Create a PyBaMM simulation set up for interation with liionpack

    Parameters
    ----------
    parameter_values : pybamm.ParameterValues class
        DESCRIPTION. The default is None.
    experiment : pybamm.Experiment class
        DESCRIPTION. The default is None.
    make_inputs : bool, optional
        Changes "Current function [A]" and "Total heat transfer coefficient
        [W.m-2.K-1]" to be inputs that are controlled by liionpack.
        The default is False.

    Returns
    -------
    sim : pybamm.Simulation
        A simulation that can be solved individually or passed into the
        liionpack solve method

    """
    # Create the pybamm model
    model = pybamm.lithium_ion.SPMe(
        options={
            "thermal": "lumped",
        }
    )
    # geometry = model.default_geometry
    if parameter_values is None:
        # load parameter values and process model and geometry
        chemistry = pybamm.parameter_sets.Chen2020
        parameter_values = pybamm.ParameterValues(chemistry=chemistry)
    # Change the current function to be an input as this is set by the external circuit
    if make_inputs:
        parameter_values.update(
            {
                "Current function [A]": _current_function,
            }
        )
        parameter_values.update(
            {
                "Current": "[input]",
                "Total heat transfer coefficient [W.m-2.K-1]": "[input]",
                "Initial concentration in electrolyte [mol.m-3]":"[input]",
            },
            check_already_exists=False,
        )

    solver = pybamm.CasadiSolver(mode="safe")
    sim = pybamm.Simulation(
        model=model,
        experiment=experiment,
        parameter_values=parameter_values,
        solver=solver,
    )
    return sim

if __name__ == "__main__":
    sim = create_simulation()
    sim.solve([0, 1800])
    sim.plot()

Changed lines 86, 122,163,256,266 in solver_utils.py.

# -*- coding: utf-8 -*-
"""
Created on Thu Sep 23 10:44:31 2021

@author: Tom
"""
import casadi
import pybamm
import numpy as np
import time as ticker
import liionpack as lp

def _mapped_step(model, solutions, inputs_dict, integrator, variables, t_eval):
    r"""
    Internal function to process the model for one timestep in a mapped way.
    Mapped versions of the integrator and variables functions should already
    have been made.

    Parameters
    ----------
    model : pybamm.Model
        The built model
    solutions : list of pybamm.Solution objects for each battery
        Used to get the last state of the system and use as x0 and z0 for the
        casadi integrator
    inputs_dict : list of inputs_dict objects for each battery
        DESCRIPTION.
    integrator : mapped casadi.integrator
        Produced by _create_casadi_objects
    variables : mapped variables evaluator
        Produced by _create_casadi_objects
    t_eval : float array of times to evaluate
        Produced by _create_casadi_objects

    Returns
    -------
    sol : list
        solutions that have been stepped forward by one timestep
    var_eval : list
        evaluated variables for final state of system

    """
    len_rhs = model.concatenated_rhs.size
    N = len(solutions)
    if solutions[0] is None:
        # First pass
        x0 = casadi.horzcat(*[model.y0[:len_rhs] for i in range(N)])
        z0 = casadi.horzcat(*[model.y0[len_rhs:] for i in range(N)])
    else:
        x0 = casadi.horzcat(*[sol.y[:len_rhs, -1] for sol in solutions])
        z0 = casadi.horzcat(*[sol.y[len_rhs:, -1] for sol in solutions])
    # t_min = [0.0]*N
    t_min = 0.0
    inputs = []
    for temp in inputs_dict:
        inputs.append(casadi.vertcat(*[x for x in temp.values()] + [t_min]))
    ninputs = len(temp.values())
    inputs = casadi.horzcat(*inputs)
    # p = casadi.horzcat(*zip(inputs, external_variables, [t_min]*N))
    # inputs_with_tmin = casadi.vertcat(inputs, np.asarray(t_min))
    # Call the integrator once, with the grid
    timer = pybamm.Timer()
    tic = timer.time()
    casadi_sol = integrator(x0=x0, z0=z0, p=inputs)
    integration_time = timer.time()
    nt = len(t_eval)
    xf = casadi_sol["xf"]
    # zf = casadi_sol["zf"]
    sol = []
    xend = []
    for i in range(N):
        start = i * nt
        y_sol = xf[:, start:start + nt]
        xend.append(y_sol[:, -1])
        # Not sure how to index into zf - need an example
        sol.append(pybamm.Solution(t_eval, y_sol, model, inputs_dict[i]))
        sol[-1].integration_time = integration_time
    toc = timer.time()
    lp.logger.debug(f"Mapped step completed in {toc - tic}")
    xend = casadi.horzcat(*xend)
    var_eval = variables(0, xend, 0, inputs[0:ninputs, :])
    return sol, var_eval

def _create_casadi_objects(I_init, htc,initial_eleConc, sim, dt, Nspm, nproc, variable_names):
    r"""
    Internal function to produce the casadi objects in their mapped form for
    parallel evaluation

    Parameters
    ----------
    I_init : float
        initial guess for current of a battery (not used for simulation).
    htc : float
        initial guess for htc of a battery (not used for simulation).
    sim : pybamm.Simulation
        A PyBaMM simulation object that contains the model, parameter_values,
        solver, solution etc.
    dt : float
        The time interval for a single timestep. Fixed throughout the simulation
    Nspm : int
        Number of individual batteries in the pack.
    nproc : int
        Number of parallel processes to map to.
    variable_names : list
        Variables to evaluate during solve. Must be a valid key in the
        model.variables

    Returns
    -------
    integrator : mapped casadi.integrator
        Solves an initial value problem (IVP) coupled to a terminal value
        problem with differential equation given as an implicit ODE coupled
        to an algebraic equation and a set of quadratures
    variables_fn : mapped variables evaluator
        evaluates the simulation and output variables. see casadi function
    t_eval : float array of times to evaluate
        times to evaluate in a single step, starting at zero for each step

    """
    inputs = {"Current": I_init, "Total heat transfer coefficient [W.m-2.K-1]": htc,"Initial concentration in electrolyte [mol.m-3]": initial_eleConc}
    solver = sim.solver
    # solve model for 1 second to initialise the circuit
    t_eval = np.linspace(0, 1, 2)
    # Initial solution - this builds the model behind the scenes
    sim.solve(t_eval, inputs=inputs)
    # step model
    # Code to create mapped integrator
    t_eval = np.linspace(0, dt, 11)
    t_eval_ndim = t_eval / sim.model.timescale.evaluate()
    inp_and_ext = inputs
    # No external variables - Temperature solved as lumped model in pybamm
    # External variables could (and should) be used if battery thermal problem
    # Includes conduction with any other circuits or neighboring batteries
    # inp_and_ext.update(external_variables)

    integrator = solver.create_integrator(
        sim.built_model, inputs=inp_and_ext, t_eval=t_eval_ndim
    )
    integrator = integrator.map(Nspm, "thread", nproc)

    # Variables function for parallel evaluation
    casadi_objs = sim.built_model.export_casadi_objects(variable_names=variable_names)
    variables = casadi_objs["variables"]
    t, x, z, p = (
        casadi_objs["t"],
        casadi_objs["x"],
        casadi_objs["z"],
        casadi_objs["inputs"],
    )
    variables_stacked = casadi.vertcat(*variables.values())
    variables_fn = casadi.Function("variables", [t, x, z, p], [variables_stacked])
    variables_fn = variables_fn.map(Nspm, "thread", nproc)
    return integrator, variables_fn, t_eval

def solve(
    netlist=None,
    parameter_values=None,
    experiment=None,
    I_init=1.0,
    htc=None,initial_eleConc=None,
    initial_soc=0.5,
    nproc=12,
    output_variables=None,
):
    r"""
    Solves a pack simulation

    Parameters
    ----------
    netlist : pandas.DataFrame
        A netlist of circuit elements with format. desc, node1, node2, value.
        Produced by liionpack.read_netlist or liionpack.setup_circuit
    parameter_values : pybamm.ParameterValues class
        A dictionary of all the model parameters
    experiment : pybamm.Experiment class
        The experiment to be simulated. experiment.period is used to
        determine the length of each timestep.
    I_init : float, optional
        Initial guess for single battery current [A]. The default is 1.0.
    htc : float array, optional
        Heat transfer coefficient array of length Nspm. The default is None.
    initial_soc : float
        The initial state of charge for every battery. The default is 0.5
    nproc : int, optional
        Number of processes to start in parallel for mapping. The default is 12.
    output_variables : list, optional
        Variables to evaluate during solve. Must be a valid key in the
        model.variables

    Raises
    ------
    Exception
        DESCRIPTION.

    Returns
    -------
    output : ndarray shape [# variable, # steps, # batteries]
        simulation output array

    """

    if netlist is None or parameter_values is None or experiment is None:
        raise Exception("Please supply a netlist, paramater_values, and experiment")

    # Get netlist indices for resistors, voltage sources, current sources
    Ri_map = netlist["desc"].str.find("Ri") > -1
    V_map = netlist["desc"].str.find("V") > -1
    I_map = netlist["desc"].str.find("I") > -1

    Nspm = np.sum(V_map)

    protocol = lp.generate_protocol_from_experiment(experiment)
    dt = experiment.period
    Nsteps = len(protocol)

    # Solve the circuit to initialise the electrochemical models
    V_node, I_batt = lp.solve_circuit(netlist)

    sim = lp.create_simulation(parameter_values, make_inputs=True)
    lp.update_init_conc(sim, SoC=initial_soc)

    v_cut_lower = parameter_values["Lower voltage cut-off [V]"]
    v_cut_higher = parameter_values["Upper voltage cut-off [V]"]

    # The simulation output variables calculated at each step for each battery
    # Must be a 0D variable i.e. battery wide volume average - or X-averaged for 1D model
    variable_names = [
        "Terminal voltage [V]",
        "Measured battery open circuit voltage [V]",
        "Local ECM resistance [Ohm]",
    ]
    if output_variables is not None:
        for out in output_variables:
            if out not in variable_names:
                variable_names.append(out)
        # variable_names = variable_names + output_variables
    Nvar = len(variable_names)
    # Storage variables for simulation data
    shm_i_app = np.zeros([Nsteps, Nspm], dtype=float)
    shm_Ri = np.zeros([Nsteps, Nspm], dtype=float)
    output = np.zeros([Nvar, Nsteps, Nspm], dtype=float)

    # Initialize currents in battery models
    shm_i_app[0, :] = I_batt * -1

    time = 0
    # step = 0
    end_time = dt * Nsteps
    step_solutions = [None] * Nspm
    V_terminal = []
    record_times = []

    integrator, variables_fn, t_eval = _create_casadi_objects(
        I_init, htc[0],initial_eleConc[0], sim, dt, Nspm, nproc, variable_names
    )

    sim_start_time = ticker.time()

    for step in range(Nsteps):
        step_solutions, var_eval = _mapped_step(
            sim.built_model,
            step_solutions,
            lp.build_inputs_dict(shm_i_app[step, :], htc,initial_eleConc),
            integrator,
            variables_fn,
            t_eval,
        )
        output[:, step, :] = var_eval

        time += dt
        # Calculate internal resistance and update netlist
        temp_v = output[0, step, :]
        temp_ocv = output[1, step, :]
        temp_Ri = np.abs(output[2, step, :])
        shm_Ri[step, :] = temp_Ri

        netlist.loc[V_map, ("value")] = temp_ocv
        netlist.loc[Ri_map, ("value")] = temp_Ri
        netlist.loc[I_map, ("value")] = protocol[step]

        # print('Stepping time', np.around(ticker.time()-tic, 2), 's')
        if np.any(temp_v < v_cut_lower):
            print("Low V limit reached")
            break
        if np.any(temp_v > v_cut_higher):
            print("High V limit reached")
            break
        # step += 1
        if time <= end_time:
            record_times.append(time)
            V_node, I_batt = lp.solve_circuit(netlist)
            V_terminal.append(V_node.max())
        if time < end_time:
            shm_i_app[step + 1, :] = I_batt[:] * -1
    all_output = {}
    all_output["Time [s]"] = np.asarray(record_times)
    all_output["Pack current [A]"] = np.asarray(protocol[: step + 1])
    all_output["Pack terminal voltage [V]"] = np.asarray(V_terminal)
    all_output["Cell current [A]"] = shm_i_app[: step + 1, :]
    for j in range(Nvar):
        all_output[variable_names[j]] = output[j, : step + 1, :]

    toc = ticker.time()
    pybamm.logger.notice(
        "Solve circuit time " + str(np.around(toc - sim_start_time, 3)) + "s"
    )
    return all_output

The example code I ran is:

import liionpack as lp
import numpy as np
import pybamm

# Generate the netlist
netlist = lp.setup_circuit(Np=16, Ns=2, Rb=1e-4, Rc=1e-2, Ri=5e-2, V=3.2, I=80.0)

output_variables = [  
    'X-averaged total heating [W.m-3]',
    'Volume-averaged cell temperature [K]',
    'X-averaged negative particle surface concentration [mol.m-3]',
    'X-averaged positive particle surface concentration [mol.m-3]',
    'X-averaged electrolyte concentration [mol.m-3]',
    ]

# Heat transfer coefficients
htc = np.ones(32) * 10
#initial_eleConc=np.ones(32) * 1000.0
initial_eleConc=np.array([1000., 1000., 2000., 2000., 1000., 1000., 1000., 1000., 1000.,
       1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000.,
       1000., 1000., 1000., 1000., 1000., 3000., 3000., 3000., 1000.,
       1000., 1000., 1000., 1000., 1000.])

# Cycling experiment, using PyBaMM
experiment = pybamm.Experiment(
    ["Charge at 50 A for 30 minutes", "Rest for 15 minutes", "Discharge at 50 A for 30 minutes", "Rest for 30 minutes"],
    period="10 seconds",
)

# PyBaMM parameters
chemistry = pybamm.parameter_sets.Chen2020
parameter_values = pybamm.ParameterValues(chemistry=chemistry)

# Solve pack
output = lp.solve(netlist=netlist,
                  parameter_values=parameter_values,
                  experiment=experiment,
                  output_variables=output_variables,
                  htc=htc,initial_eleConc=initial_eleConc)
lp.plot_output(output)

However, in the simulation the first value of initial electrolyte concentration, i.e., initial_eleConc[0] is populating across all the cells. Here is the output showing that : I think the issue is in line 256 of solver_utils.py, i.e.

integrator, variables_fn, t_eval = _create_casadi_objects(
        I_init, htc[0],initial_eleConc[0], sim, dt, Nspm, nproc, variable_names
    )

where it is only taking the first element of the array initial_eleConc[0]. I gave initial_eleConc array similar to htc(heat transfer coefficient) array. As htc array is going into lumped model, so every cell can have their own individual heat transfer coefficient. In contrast, the parameters in the base model seems not populating across the cells.

[valentinsulzer]

Haven't looked specifically into what's happening here, but we should come up with a user interface that allows users to change any parameter like this, instead of adding each one as another argument to all the functions

[TomTranter]

Can you make a branch from main with the code changes you made?

[TomTranter]

Haven't looked specifically into what's happening here, but we should come up with a user interface that allows users to change any parameter like this, instead of adding each one as another argument to all the functions

Agreed. The current has to be an input for Liionpack to work and is changed on every step. HTC technically doesn't but you would want it to be if solving a dynamic thermal problem and again might change a little bit for every step. The rest are more like ensemble constants and would be different to start with but then not change. Is it worth making this distinction and calling it something else?

[ksnvikrant]

Can you make a branch from main with the code changes you made?

I created a branch name issue-57-init_cond with the above changes

[TomTranter]

Hi @ksnvikrant this issue is also related to being able to pass in the model flexibly. I think the answer to both is to set up the simulation and then pass that into the liionpack solvers. Only problem there is how to enforce the current to be an input. I guess throwing errors if it's not might be ok for now. The code changes you made will now be out of date after the solver refactor so I think it's probably easier just to make them again off a new branch. I can combine this with the DFN model issue into a new branch

[TomTranter]

@ksnvikrant I've just pushed a bunch of changes which should help you with this issue. There is a new notebook on using inputs. Let me know if you have any questions

[ksnvikrant]

Hi @TomTranter Thanks for the changes. I tried to create an example for different initial SOC condition for the cells in pack. However, it is still taking the first initial SOC of the cell and populating across all other cells. In order to run the example below, first you need to comment line 44 of solvers.py, i.e., #lp.update_init_conc(self.simulation, SoC=initial_soc). The example code is here:

import liionpack as lp
import pybamm
import numpy as np
I_mag = 5.0
OCV_init = 4.0  # used for intial guess
Ri_init = 5e-2  # used for intial guess
R_busbar = 1.5e-6  # very small to simulate evenly distributed currents
R_connection = 1e-2
Np = 4
Ns = 1
Nbatt = Np * Ns
netlist = lp.setup_circuit(
Np=Np, Ns=Ns, Rb=R_busbar, Rc=R_connection, Ri=Ri_init, V=OCV_init, I=I_mag
)

experiment = pybamm.Experiment( ["Discharge at 5 A for 30 minutes", "Rest for 15 minutes"], period="10 seconds", )

output_variables = [
'X-averaged negative particle surface concentration [mol.m-3]', 'X-averaged positive particle surface concentration [mol.m-3]', ]

chemistry = pybamm.parameter_sets.Chen2020 parameter_values = pybamm.ParameterValues(chemistry=chemistry) socInit=np.array([0.2, 0.3, 0.4, 0.5])

c_n_max = parameter_values["Maximum concentration in negative electrode [mol.m-3]"] c_p_max = parameter_values["Maximum concentration in positive electrode [mol.m-3]"]

neg_init=np.zeros(NpNs) pos_init=np.zeros(NpNs)

for i in range(NpNs): x, y = pybamm.lithium_ion.get_initial_stoichiometries(socInit[i], parameter_values) neg_init[i], pos_init[i] = x c_n_max, y * c_p_max

parameter_values.update({"Initial concentration in negative electrode [mol.m-3]": "[input]"}) parameter_values.update({"Initial concentration in positive electrode [mol.m-3]": "[input]"}) inputs = {"Initial concentration in negative electrode [mol.m-3]": neg_init,"Initial concentration in positive electrode [mol.m-3]": pos_init}

print (np.unique(neg_init))

output = lp.solve( netlist=netlist, parameter_values=parameter_values, experiment=experiment, inputs=inputs,output_variables=output_variables ) lp.plot_output(output) lp.show_plots()

The output for surface concentration of negative and positive particles for the cells should look different as the initial SOC of cells are different. However, output is same for all the cells. 
![Figure_5](https://user-images.githubusercontent.com/92038470/143487403-912c4bda-e3d0-4f06-96f8-d7f6eb76b246.png)
![Figure_6](https://user-images.githubusercontent.com/92038470/143487426-a3e34ef9-4133-4a59-acea-b0f72e361ea5.png)

[TomTranter]

We could change the logic in the solvers to not use the initial_soc argument if is None and change the default and that should allow you to customize in this way or we could develop the logic to accept an array of initial_soc

[TomTranter]

This issue should now be addressed by PR #103 which is merged. See https://github.com/pybamm-team/liionpack/blob/main/examples/different_initial_soc.py for example and let me know if we can close this issue

[TomTranter]

@ksnvikrant also that other branch is redundant now so can you delete it?