Can I get a better fit to data

I've written a code to find some parameters to a mathematical model but I can't seem to get a good fit to the experimental data. I've seen that this model can model the data very well but I can't seem to get the right parameters that actually give me a reasonable fit. Are there any tricks I can use to get the right parameters to get me close to the data?

The functions and code I use are:

function y=eta_SPM(c_s,k,T,I_app)
if isrow(I_app)==true
    I_app=I_app';
end
if isrow(c_s)==true
    c_s=c_s';
end
R=8.314;
F=9.648e4;
g_0=2*k*sqrt(c_s.*(1-c_s));
y=(2*R*T/F)*asinh(I_app./g_0);

.

function V=volt(X,L_a,L_c,T,I_app,t)
%Compute the lithium concentration in each of the 
r=linspace(0,1,100);
D_a=X(1);
D_c=X(2);
u_0_a=X(3);
u_0_c=X(4);
k_a=X(5);
k_c=X(6);
alpha_a=X(7);
alpha_c=X(8);
c_c=c(u_0_c,alpha_c*I_app,D_c,r,t);
c_a=c(u_0_a,-alpha_a*I_app,D_a,r,t);

%Compute Voltage
V=OCV_plus(c_c)+eta_SPM(c_c,k_c,T,-alpha_c*I_app/L_c)-(OCV_minus(c_a)+eta_SPM(c_a,k_a,T,alpha_a*I_app/L_a));

.

function y=OCV_plus(x)
if isrow(x)==true
    x=x';
end
y=-2.613*x.^6+9.858*x.^5-16.63*x.^4+14.09*x.^3-4.952*x.^2-0.4427*x+4.27;

.

function y=OCV_minus(x)
if isrow(x)==true
    x=x';
end
y=(289.7*x.^2+339.3*x+99.95)./(x.^3+4408*x.^2+4080*x+142.9);

Here I use a weight function to highlight the parts of the data I want the optimiser to take special notice of.

function y=objective_fn(V_exp,V_calc,t)

%This is the objective function
N=length(t); %gets length of vectors
if isrow(V_exp)==true
    V_exp=V_exp';
end
if isrow(V_calc)==true
    V_calc=V_calc';
end

if isrow(t)==true
    t=t';
end
%define the weight of the objective function to concentrate on a particular
%feature
w=ones(N,1);
w(6500:6800)=2;
w(1:100)=2;
y=trapz(t,(V_exp.*w-V_calc).^2);

.

%This program finds the optimised parameters 

%Define upper and lower bounds for the parameters;
D_a=0.03;
D_c=0.01;
u_0_a=0.95;
u_0_c=0.27;
k_a=1e-3;
k_c=1e2;
alpha_a=0.01;
alpha_c=0.015;
lb = [10^-4, 10^-4 ,0.7 ,0.15 ,1e-3 ,1e-3 ,1e-4 ,1e-4 ];
ub = [1e1, 1e1, 0.99, 0.5, 2 ,2 ,1 ,1 ];
x0=[D_a,D_c,u_0_a,u_0_c,k_a,k_c,alpha_a,alpha_c];
n=length(t);
T=293;
I_app=ones(1,n);
L_a=1;
L_c=1;
V_0=volt(x0,L_a,L_c,T,I_app,t);
A = [];
b = [];
Aeq = [];
beq = [];
fun=@(X) objective_fn(V_exp,  volt(X,L_a,L_c,T,I_app,t)  ,t);
X = fmincon(fun,x0,A,b,Aeq,beq,lb,ub);
%Compare with experimental data.

D_a=X(1);
D_c=X(2);
u_0_a=X(3);
u_0_c=X(4);
k_a=X(5);
k_c=X(6);
alpha_a=X(7);
alpha_c=X(8);

r=linspace(0,1,150);
c_c=c(u_0_c,alpha_c*I_app,D_c,r,t);
c_a=c(u_0_a,-alpha_a*I_app,D_a,r,t);

%Compute Voltage
V_op=OCV_plus(c_c)+eta_SPM(c_c,k_c,T,-alpha_c*I_app/L_c)-(OCV_minus(c_a)+eta_SPM(c_a,k_a,T,alpha_a*I_app/L_a));
plot(t,V_exp);
hold on
plot(t,V_op);
xlabel('Time(hrs)');
ylabel('Terminal Voltage');
axis([0 1.2 2.7 4.5]);

I did a curve fit to the experimental data which is valid for 0<=t<=2.05

function y=V_data(t)
if isrow(t)==true
    t=t';
end
y=(-3821*t.^4+1209*t.^3-2564*t.^2+3545*t+2284)./(t.^5-657.6*t.^4*t.^3+48.25*t.^2+914.6*t+544.6);

.

function u=c(v_0,q,D,r,t)
n=length(r);m=length(t);  %The dimensions of u are defined by the length os r and t.
v=zeros(m,n); %Initialise u
dr=r(2); %as both t and r start at zero, dr and dt can be defined in this way
dt=t(2);
alpha=D*dt/dr^2; %compute alpha and beta
beta=D*dt/(2*dr);
A_main=-(1+alpha)*ones(n,1); %Begin constructing the vectors to put in the diagonals
A_u=0.5*alpha+beta./r(1:n-1);
A_l=0.5*alpha-beta./r(2:n);
A=diag(A_u,1)+diag(A_l,-1)+diag(A_main,0); %Construct A
B_main=(alpha-1)*ones(n,1); %Only the main diagonal for differes by anyting more than a sign
B=diag(-A_u,1)+diag(B_main,0)+diag(-A_l,-1); %Construct B
A(1,1)=-(1+3*alpha); %Insert the inner boundary values for A and B
A(1,2)=3*alpha;
B(1,1)=3*alpha-1;
B(1,2)=-3*alpha;
A(n,n)=-(1+alpha); %Insert the outer boundary values for A and B
A(n,n-1)=alpha;
B(n,n)=-alpha;
B(n,n-1)=alpha-1;
p=zeros(n,1); %Used in the consutruction of c
p(n)=1;
B_1=A\B;
w=A\p;
gamma=-2*(dr/D)*(0.5*alpha+beta/r(n)); %For computational efficency.
v(1,:)=v_0*ones(1,n);
for i=2:m
c=gamma*w*(q(i-1)+q(i));
sol=B_1*v(i-1,:)'+c;
v(i,:)=sol';
end
u=v(:,end);

[Torsten]

I've seen that this model can model the data very well but I can't seem to get the right parameters that actually give me a reasonable fit.

Is this not a contradiction ?

[Steven Lord]

I did a curve fit to the experimental data which is valid for 0<=t<=2.05

Are you sure about that? Check your denominator -- if I change t.^4*t.^3 to t.^4+t.^3 (which I suspect is what you intended to write) I see that it has a root in that range.

denominator = @(t) t.^5-657.6*t.^4+t.^3+48.25*t.^2+914.6*t+544.6;
fplot(denominator, [0 2.05])
yline(0)
r = fzero(denominator, 1.5);
xline(r, 'r')
title("Root at " + r)

Let's check the value of the numerator at that point.

numerator = @(t) -3821*t.^4+1209*t.^3-2564*t.^2+3545*t+2284;
numerator(r)

ans = -5.3266e+03

numerator(r)./denominator(r)

ans = -4.6853e+16

Does your data have a -Inf (or a large magnitude negative value) at that point?

[Matthew Hunt]

Hi, I see where the issue is and the proper code for the curve fit of the experimental data is:

function y=V_data(t)
if isrow(t)==true
    t=t';
end
y=(-2.376e+04*t.^4+5.47e+04*t.^3-3.137e+04*t.^2+3.294e+04*t+2.165e+04)./(t.^5-6368*t.^4+1.418e+04*t.^3-6837*t.^2+8408*t+5154);

The input is in hours and it goes from 0 to 2.03.

Answer 1

[Yoga]

Hi @Matthew Hunt,

I guess you found the fix for the question.

The proper code for the curve fit of the data is working fine as you have mentioned in the comments.

function y=V_data(t)
if isrow(t)==true
    t=t';
end
y=(-2.376e+04*t.^4+5.47e+04*t.^3-3.137e+04*t.^2+3.294e+04*t+2.165e+04)./(t.^5-6368*t.^4+1.418e+04*t.^3-6837*t.^2+8408*t+5154);

[Matthew Hunt]

Actually I didn't.