Can I get a better fit to data
3 views (last 30 days)
Show older comments
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);
3 Comments
Steven Lord
on 2 Feb 2023
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)
numerator(r)./denominator(r)
Does your data have a -Inf (or a large magnitude negative value) at that point?
Answers (1)
Yoga
on 10 Mar 2023
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);
See Also
Categories
Find more on Nonlinear Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!