Write the transfer function of the circuit that models the complex impedance of the cell.

tf=(Rser+(Rfast/((Rfast s Cfast)+1))+(Rslow/((Rslow s Cslow)+1)))

Convert the tf into a transfer function model.

tfm=TransferFunctionModel[{{tf}},s]

Generate a battery current profile.

Ibatt=UnitStep[t-1]-UnitStep[t-5]

Plot[Ibatt,{t,0,10},Exclusions->None]

Assign values to the circuit parameters.

vals={Rser->0.1,Rfast->0.05,Cfast->10,Rslow->0.01,Cslow->100}

Voc=4.2

Find the output response for these circuit values and current profile.

y=OutputResponse[tfm/.vals,Ibatt,{t,10}];

Plot[Voc-y,{t,0,10}]

Define some different circuit values to get a different response and plot both responses.

dvals={Rser->0.06,Rfast->0.08,Cfast->15,Rslow->0.015,Cslow->50}

dy=OutputResponse[tfm/.dvals,Ibatt,{t,10}];

Plot[{Voc-y,Voc-dy},{t,0,10}]

Find the area of the difference between the two responses.

Plot[Abs[y-dy],{t,0,10}]

NIntegrate[Abs[y-dy],{t,0,10}]

Given only the first response curve as a template, find the circuit values that produced that response as a test of the modeling algorithm. Do this by minimizing the difference between the template and the output of a model whose circuit values are adjusted by NMinimize.

First write a dummy function which will force NMinimize to use numeric values for the objective function parameters instead of trying to evaluate the function symbolically.

f[{Rse_?NumericQ,Rfa_?NumericQ,Cfa_?NumericQ,Rsl_?NumericQ,Csl_?NumericQ}]:=(error=NIntegrate[Abs[OutputResponse[tfm/.{Rser->Rse,Rfast->Rfa,Cfast->Cfa,Rslow->Rsl,Cslow->Csl},Ibatt,{t,0,10}]-y
],{t,0,10}
];
First[error]
);

f[{0.11,0.051,11,0.011,110}]

Now fit the model and time how long it takes.

sol=Timing[NMinimize[{f[{Rser,Rfast,Cfast,Rslow,Cslow}],0<Rser<0.5,0<Rfast<0.5,10<Cfast<100,0<Rslow<0.5,50<Cslow<200},{Rser,Rfast,Cfast,Rslow,Cslow},Method->"DifferentialEvolution"]]

sol[[2,2]]

Plot the optimized solution along with the known correct solution.

yopt=OutputResponse[tfm/.sol[[2,2]],Ibatt,{t,10}];

Plot[{Voc-y,Voc-yopt},{t,0,10}]

Plot[Abs[y-yopt],{t,0,10},PlotRange->All]

The open-circuit voltage sags during discharge, so we have to include that in the model.

sag=Piecewise[{{0,t<1},{(0.05/4) (t-1),1≤t≤5},{0.05,t>5}}]

Plot[sag,{t,0,10}]

Plot[(Voc-sag)-y,{t,0,10}]

Now attempt to fit a model to real measurement data. First import the data.

data3p3=Import["/Users/mcgill/Documents/My\ Projects/Lithium\ Battery\ Safety/Test\ Data/NCR18650B/3300mA/NCR18650B_3300mA.xls"];

Vbat=Drop[data3p3[[1,All,3]],1];

Length[Vbat]

ListPlot[Vbat,PlotRange->{{0,12500},{2.5,4.5}},GridLines->{Range[0,12000,1000],Range[2.5,4.5,0.5]},Frame->True,FrameStyle->Directive[18],FrameLabel->{Style["Time (secs)",FontSize->18],Style["Cell Voltage (V)",FontSize->18]},PlotLabel->Style["NCR18650B Discharge\r6 min at 3.3 A, 15 min rest",FontSize->24],PlotStyle->Red]

ListPlot[Vbat[[800;;2100]],PlotRange->{{0,1000},{3.7,4.2}},GridLines->{Range[0,1300,200],Range[3.7,4.2,0.1]},Frame->True,FrameStyle->Directive[18],FrameLabel->{Style["Time (secs)",FontSize->18],Style["Cell Voltage (V)",FontSize->18]},PlotLabel->Style["NCR18650B Discharge\r6 min at 3.3 A, 15 min rest",FontSize->24],PlotStyle->Red]

Vbat[[890;;910]]

Vbat[[1250;;1270]]

Vbatp1=Vbat[[800;;1799]];

Length[Vbatp1]

Vbatp1[[95;;105]]

Vbatp1[[455;;465]]

Ibatp1=Piecewise[{{0,t<100},{3.3,100≤t≤458},{0,t>458}}];

Plot[Ibatp1,{t,0,999}]

{Vbatp1[[1]],Vbatp1[[999]]}

Vocp1=Table[Piecewise[{{Vbatp1[[1]],t<100},{Vbatp1[[1]]-((Vbatp1[[1]]-Vbatp1[[999]])/359) (t-100),100≤t≤458},{Vbatp1[[999]],t>458}}],{t,0,999}];

ListPlot[Vocp1]

ListPlot[Flatten[Table[Voc-y,{t,0,10,0.01}]]]

ManhattanDistance[Flatten[Table[y,{t,0,10,0.01}]],Flatten[Table[dy,{t,0,10,0.01}]]]

Timing[ManhattanDistance[Flatten[Table[y,{t,0,10,0.01}]],Flatten[Table[dy,{t,0,10,0.01}]]]]

g[{Rse_?NumericQ,Rfa_?NumericQ,Cfa_?NumericQ,Rsl_?NumericQ,Csl_?NumericQ}]:=Module[{modResp,error},
modResp=OutputResponse[tfm/.{Rser->Rse,Rfast->Rfa,Cfast->Cfa,Rslow->Rsl,Cslow->Csl},Ibatp1,{t,0,999}];
error=ManhattanDistance[Flatten[Table[modResp,{t,1000}]],Vbatp1];
First[error];
]; (* end module *)

bigvals={Rser->0.05,Rfast->0.05,Cfast->500,Rslow->0.05,Cslow->8000}

modResp=OutputResponse[tfm/.bigvals,Ibatp1,{t,0,999}];

ListPlot[Vocp1-Flatten[Table[modResp,{t,0,999}]]]

ListPlot[Vocp1-Flatten[Table[OutputResponse[tfm/.{Rser->0.05,Rfast->0.05,Cfast->500,Rslow->0.05,Cslow->8000},Ibatp1,{t,0,999}],{t,0,999}]]]

test=Flatten[Table[modResp,{t,0,999}]];

Length[test]

Vsagp1=Table[Piecewise[{{Vbatp1[[1]],t<100},{Vbatp1[[1]]-((Vbatp1[[1]]-Vbatp1[[999]])/359) (t-100),100≤t≤458},{Vbatp1[[999]],t>458}}],{t,0,999}];

test[[1;;10]]

test2=test-Vsagp1;

Length[test2]

test2[[1;;10]]

ListPlot[Vsagp1]

Ibatp1k=Table[Piecewise[{{0,t<100},{3.3,100≤t≤458},{0,t>458}}],{t,0,999}];

ListPlot[Ibatp1k]

dtm=ToDiscreteTimeModel[tfm,1,z]

outk=Flatten[OutputResponse[dtm/.bigvals,Ibatp1k]];

Length[outk]

ListPlot[Vsagp1-outk]

Manipulate[ListPlot[Vsagp1-Flatten[OutputResponse[dtm/.{Rser->a,Rfast->b,Cfast->c,Rslow->d,Cslow->e},Ibatp1k]
],PlotRange->{{0,1000},{3.0,4.2}},GridLines->Automatic
],{{a,0.05,"Rser"},0,0.1},{{b,0.05,"Rfast"},0,0.1},{{c,1000,"Cfast"},0,2000},{{d,0.05,"Rslow"},0,0.1},{{e,10000,"Cslow"},0,20000}]