Can I get a better fit to data

2 次查看（过去 30 天）
Matthew Hunt
Matthew Hunt 2023-2-2
评论： Matthew Hunt 2023-3-14
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 个评论
显示 1更早的评论
Steven Lord
Steven Lord 2023-2-2
Ran in:
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
Matthew Hunt 2023-2-3
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Yoga
Yoga 2023-3-10
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);

类别

在 Help Center 和 File Exchange 中查找有关 Data Clustering 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by