lsqcurvefit doesn't curve fit

I have a model which I want to paramtrise using lsqcurvefit. I have 10 parameters that I must find and I have 10 pieces of data (or more) that I can call on. I set up my function that I want to minimise including the function which includes the model. When I use 10 points I get the message that a minimum is possible and when I plot the solution using the parameters and compare it against the experimental data, I get completely different curves, the solution should overly the points I get but that just isn't the case.

Any idea why this would happen?

Mat

Answer 1

[Star Strider]

Use as many data as you have. Also, nonlinear parameter estimation techniques are very sensitive to the initial estimates (that you give to the routine to start with), and an inaccurate set can cause the routine to end up in a local minimum rather than a minimum that is much closer to the correct parameters. Choosing the correct values can be challenging.

If you repeatedly have problems guessing the correct initial parameter values, use one of the Global Optimization Toolbox functions (such as the genetic algorithm ga function) to search out the best parameter set. Those take time, however they are usually succesful. (For ga, begin with a large initial population, so it has a better probability of discovering the best parameter set.)

[Matthew Hunt]

I'm loathed to use ga algorithms as they can take forever to get something. I hear you about a local minimum.

My first idea was to treat the problem as a nonlinear system of algebraic equations, this should be rather quick to do I can use Newton Raphson. That didn't work particularly well, and I'm not too sure why.

Mat

[Star Strider]

The MATLAB ga alfgorithm can actually be reasonably fast, although with 10 parameters, it’s best you find something else to do while it’s running.

Newton-Raphson is a root-finding algorithm. If you’re estimating a nonlinear system of algebraic equations and you want to find a set of parameters that will equate your equations to 0, use the fsolve function.

I’m still not sure what you’re doing. Expanding on a description of that could make it easier to converge on a solution.

[Matthew Hunt]

I want to do something more clever that just throw a genetic algorithm at the problem.

I think I can at least get a good starting point for the other algorithms by treating the parametrisation problem as a nonlinear algebraic equation problem. The model I want to fit, is quite a simpe one as it's analytically solveable. My first thought was to evaluat the model at the same points as the dataset and set this up as a nonlinear algebraic system of equations.

After that, possibly use the results as an initial guess for a more sophisticated method.

[Matt J]

My first thought was to evaluat the model at the same points as the dataset and set this up as a nonlinear algebraic system of equations.

There is little or no distinction between that and what lsqcurvefit (or fsolve) is already doing. And it doesn't solve the problem because nonlinear equation solvers also rely on accurate initial parameter guesses, as you probably discovered when you tried Newton-Raphson.

It would help a lot if we could see the actual model function and data.

[Matthew Hunt]

Unforttunately I would rather show my data.

[Star Strider]

@Matthew Hunt —

‘I want to do something more clever that just throw a genetic algorithm at the problem.’

A genetic algorithm is clever enough, if you have no idea what the correct initial parameter set should be.

Meanwhile, you still haven’t enlightened us as to what your actual problem is. We can keep guessing at a solution hoping (against all odds) to come up with the correct one by sheer chance, or you can let us help you with it.

[Matt J]

I want to do something more clever that just throw a genetic algorithm at the problem...The model I want to fit, is quite a simpe one as it's analytically solveable.

@Matthew Hunt,

If the fit has an analytical solution, then using lsqcurvefit, or any other iterative solver, is not clever. You should just code the analytical solution algorithm.

[Matthew Hunt]

I'm trying to fit parameters to the single particle model in lithium ion battery modelling. It involves solving a spherical diffusion equation (which has an analytical solution). I'm currently using the terminal velocity as experimental data, so I'm using expressions for overpotential.

[Matthew Hunt]

The model I am using is here:

https://arxiv.org/pdf/1702.02471.pdf

[Star Strider]

Thanks for the paper. It’s in my clip file.

It would be easier to work with your code, since that seems to be where the problem is, and I can’t figure out what expression in that paper you’re using for your model.

Post your code, describe what it’s supposed to do in some detail, and describe what you want to do with it. Is it an optimization problem without data, a parameter estimation problem with data, or something else?

[Matthew Hunt]

%This is the code for the inverse problem for the single particle model.

%The inverse problem will use the terminal voltage for the inverse problem

%This computes the constants for the Greens function:
N=15; %Number of terms used in Green's function
d=0.01; %Using the small different in the intevals to avoid infinity
f=@(x)tan(x)-x; %The transecendal equation 
mu_n=zeros(N,1); %array which stores the solutions
mu_n(1)=fzero(f,pi/2-d);
for i=0:N-1
    mu_n(i+2)=fzero(f,[(2*i+1)*pi/2+d (2*(i+1)+1)*pi/2-d]);
end
mu_n=mu_n(2:N+1);

%Import data as time and voltage:
load VoltageData.mat;
load('OCP_Cathode.mat'); %Load the OCP data for the cathode
load('OCP_aNODE.mat'); %Load OCP data for the anode;
SOC_a=OCP_Anode(:,1);
SOC_c=OCP_Cathode(:,1);
OCV_a=flip(OCP_Anode(:,2));
OCV_c=OCP_Cathode(:,2);
time=VoltageData(:,1);
voltage=VoltageData(:,2); 
%Import the applied current data
%Set it to the correct value:
I_app=ones(1,10);
Time=time/time(end);
%Get the data to be used:
t_data=linspace(0,Time(end),10);
L=length(Time);
V_data=interp1(Time,voltage,t_data);
%Choose initial conditions for the parameters:
X_0=0.5*ones(1,10);
X_0(5:10)=10^-7;

lb=zeros(1,10);
ub=ones(1,10);
%Set up the function for lsqcurvefit:

fun = @(X,t)terminal_voltage(V_data,I_app,mu_n,t_data,X,SOC_a,SOC_c,OCV_c,OCV_a)';
options=optimset('MaxFunEvals',100000,'TolX',10^-10);
x = lsqcurvefit(fun,X_0,t_data',V_data,lb,ub,options);
X=x;

%Compare the theoretical prediction to experimental data
SOC_anode=anode_c(I_app,mu_n,1,t_data,X);
SOC_cathode=cathode_c(I_app,mu_n,1,t_data,X);

O_plus=interp1(1-SOC_c,OCV_c,SOC_cathode);
O_minus=interp1(SOC_a,OCV_a,SOC_anode);
V=(eta_plus(SOC_cathode,I_app,X)-eta_minus(SOC_anode,I_app,X)+O_plus-O_minus);
plot(t_data,V);
hold on
plot(t_data,V_data,'ro');
title('Comparison of Theory and Experiment');
legend('Calculated Voltage', 'Experimental Voltage');
xlabel('Time');
ylabel('Terminal Voltage');

[Matthew Hunt]

function C = anode_c(I,mu_n,r,t,X)
C_1=zeros(length(t),1);
C_2=zeros(length(t),1);
c_0=X(2);
D=X(4);
Gamma=X(6);
gamma=X(8);
nu=X(10);
I_bar=max(abs(I));
I_hat=I/I_bar;
%This is the solution of the spherical diffusion equation using Greens
%Function.
if (length(c_0)==1)
    N_r=500; %Steps used in the r_bar integration
else
    N_r=length(c_0);
end
N_t=length(I); %Steps used in the tau integration.
r_bar=linspace(0,1,N_r);
G_1=Greens_fn(D,mu_n,r,t,r_bar);
for i=2:N_t
    C_1(i)=trapz(r_bar,c_0.*G_1(i,:));
    G_2=Greens_fn(D,mu_n,1,t(i)-t(1:i),1);
    C_2(i)=trapz(t(1:i),I_hat(1:i)'.*G_2);
end
C=C_1+(D*Gamma)*C_2;
C(1)=c_0;
end

[Matthew Hunt]

function C = cathode_c(I,mu_n,r,t,X)
C_1=zeros(length(t),1);
C_2=zeros(length(t),1);
C_test=zeros(length(t),1);
c_0=X(1);
D=X(3);
Gamma=X(5);
gamma=X(7);
nu=X(9);
I_bar=max(abs(I));
I_hat=I/I_bar;
%This is the solution of the spherical diffusion equation using Greens
%Function.
if (length(c_0)==1)
    N_r=500; %Steps used in the r_bar integration
else
    N_r=length(c_0);
end
N_t=length(I); %Steps used in the tau integration.
r_bar=linspace(0,1,N_r);
G_1=Greens_fn(D,mu_n,r,t,r_bar);
for i=2:N_t
    C_1(i)=trapz(r_bar,c_0.*G_1(i,:));
    G_2=Greens_fn(D,mu_n,1,t(i)-t(1:i),1);
    C_2(i)=trapz(t(1:i),I_hat(1:i)'.*G_2);
end
C=C_1-(D*Gamma)*C_2;
C(1)=c_0;
end

[Matthew Hunt]

function e_m = eta_minus(c_anode,I,gamma_minus,nu_minus)
%This is the overpotential function.
i_0=sqrt(c_anode.*(1-c_anode));
e_m=asinh(gamma_minus*I'./(nu_minus*i_0));
end

function e_p = eta_plus(c_cathode,I,gamma_plus,nu_plus)
%This is the overpotential function.
i_0=sqrt(c_cathode.*(1-c_cathode));
e_p=asinh(gamma_plus*I'./(nu_plus*i_0));
end

[Matthew Hunt]

function G = Greens_fn(D,mu_n,r,t,r_bar)
%This is the Green's function for spherical diffusion equation
%Solution can be found in A.D. Polyanin. The Handbook of Linear Partial Differential Equations
%for Engineers and Scientists, Chapman & Hall, CRC 17
G=zeros(length(t),length(r_bar));
N_t=length(t);
N=length(mu_n); %Number of terms used in Green's function series
for i=1:N
    for j=1:N_t
    G_i(j,:)=(2*r_bar/r).*((1+mu_n(i)^-2).*sin(mu_n(i)*r).*sin(mu_n(i)*r_bar).*exp(-D*mu_n(i)^2*t(j))); %Ther series part of the Greens function
    G(j,:)=G(j,:)+G_i(j,:); %Adding uo the terms of the series
    end
end
for i=1:N_t
    G(i,:)=G(i,:)+3*r_bar.^2; %Adding the final part to give the Green's function
end
end

[Matthew Hunt]

function V=terminal_voltage(V_exp,I_app,mu_n,t,X,SOC_a,SOC_c,OCV_c,OCV_a)
gamma_plus=X(5);
gamma_minus=X(6);
nu_plus=X(9);
nu_minus=X(10);
SOC_cathode=cathode_c(I_app,mu_n,1,t,X);
SOC_anode=anode_c(I_app,mu_n,1,t,X);
O_plus=interp1(1-SOC_c,OCV_c,SOC_cathode);
O_minus=interp1(SOC_a,OCV_a,SOC_anode);
V=eta_plus(SOC_cathode,I_app,gamma_plus,nu_plus)-eta_minus(SOC_anode,I_app,gamma_minus,nu_minus)+O_plus-O_minus-V_exp';

Answer 2

[Matt J]

You could have a bug in your model function, such that it is not implementing the curve you that you think it is. What happened when you used your model function code to generate a curve with known parameters? Did the curve look as expected?

[Matthew Hunt]

The model function is fine.