Optimization Curve Fitting: Curves are not converging

I'm trying to fit some experimental data with a mathemtical model. I'm trying to using Optimization techiques to do the fitting, namely "fminsearch, and fmincon". However, and as you can see in the attcahed plot image, the cuvres do not match. They do some how have simlar shaeps, and tried to follow MATLAB's guide on the topic as well.

My fitting parametres are:

• Q. (Highest Priority)
• Total (2nd Highest)
• S1 and 2.
• Tep.

I have attached my sample data along with my code. It takes a little while to run but it works. Looking forward for your help, and thank you in advance.

[Matt J]

and as you can see in the attcahed plot image

There is no image attachment at the time I'm writing this.

[Trung]

Thank you for your comment and sorry for mssing the fig. It is attached now.

[Matt J]

It takes a little while to run but it works.

The time consuming part is the calculation of Pre. I urge you to make life easy on us and just attach that, along with the other fixed problem parameters, in a .mat file. Then we can jump straight to the optimization.

[Trung]

I expremely apoligize for this. I tried to keep these two line incase you need fast calculation.

Real_Time = downsample(Experiment(:,1),20); 
Real_Pre = downsample(Experiment(:,2),20); 

[Matt J]

Thanks you for that part, but the calculation of Pre is still quite slow and it is unnecessary for your question to have us repeat it. If you have already done the part below, then just attach all these pre-calculated variables in a .mat fil (one file, not several).

Experiment = importdata('Data.txt');

Real_Time = Experiment(:,1); 
Real_Pre = Experiment(:,2); 
%If it is taking so long, this will make it faster 
% Real_Time = downsample(Experiment(:,1),20); %s
% Real_Pre = downsample(Experiment(:,2),20); %Torr

 Qp = 10; 
Total=8e16; 
Tim= Real_Time; 
S1= 15; 
S2=50; 
Tep= 293; 


Pre=Fitt_Test(Qp,Total,Tim,S1,S2,Tep);

[Trung]

Again. I really apoligize for the trouble. I have attached the run resutls as a Mat file.

Answer 1

[Matt J]

One immediate problem that I see is that in sseval, the optimization variables x(i) aren't used at all in the calculation of sse. That's what the squiggly red underlines are trying to warn you of:

[Trung]

Thank you for your comment! I have modified it as shown below:

function sse = sseval(Pre,Real_Pre)
%  Qp = x(1); 
% Total=x(2); 
% S1= x(3); 
% S2=x(4); 
% Tep= x(5);
% Ton = Tim;
% Pror = Fitt_Test(x(1),x(2),Ton,x(3),x(4),x(5));
 sse = sum((Pre - Real_Pre).^2);
% sse =  sqrt(mean((Pre-Pror).^2));
end

since the x has no use. This didnot change the results.

[Matt J]

since the x has no use

The x is the important part. It's what the solver is trying to find. So, the sse has to depend on it.

[Trung]

Thank you for your replay. When I followed the tutorial, they used x to find the see:

function sse = sseval(x,tdata,ydata)
A = x(1);
lambda = x(2);
sse = sum((ydata - A*exp(-lambda*tdata)).^2);

In my case, I thought I only need the real data(Real_Pre), and the model one (Pre).

[Matt J]

Yes but Pre must be computed inside sseval() and it must be computed based on x, which is what fmincon is passing to sseval through fun(). In your current implementation, Pre is just some fixed variable that gets passed to sseval every time fmincon calls it. We can see this by calling fun() several times with random inputs. The output is always the same:

load Fit.mat

Warning: Could not find appropriate function on path loading function handle C:\Users\aiman\OneDrive\Diffusion_Research\Outgassing\MAT\Fitting.m>@(x)sseval(Real_Pre,Pre)

fun = @(x)sseval(x,Tim,Real_Pre,Pre);

fun(rand(5,1))

ans = 7.0506e-04

fun(rand(5,1))

ans = 7.0506e-04

fun(rand(5,1))

ans = 7.0506e-04

function sse = sseval(x,Tim,Pre,Real_Pre)
Eng = x(1);
Molec = x(2);

Sped1 = x(3);
Sped2 = x(4);
Tempre = x(5);
% sse = sum((Pre - Real_Pre).^2);
sse =  sqrt(mean((Pre-Real_Pre).^2));
end

[Trung]

Thank you for you valuble input! I have implemented your remarks as shown in the code below:

Experiment = importdata('Data.txt');

% Real_Time = Experiment(:,1); 
% Real_Pre = Experiment(:,2); 
%If it is taking so long, this will make it faster 
Real_Time = downsample(Experiment(:,1),20); %s
Real_Pre = downsample(Experiment(:,2),20); %Torr

%  Qp = 10; 
% Total=9e20; 
% Tim= Real_Time; 
% S1= 15; 
% S2=50; 
% Tep= 293; 
% Pre=Fitt_Test(Qp,Total,Tim,S1,S2,Tep);

fun = @(x)sseval(x,Real_Time,Real_Pre);
x0 = [10;9e20;15;20;293];
options = optimset('TolX',1e-9,'Algorithm','sqp');%'Algorithm','sqp'
% hybridopts = optimoptions('fmincon','OptimalityTolerance',1e-1);
% options = optimoptions('ga','HybridFcn',{'fmincon',hybridopts});
lb = [3,7e16,3,30,290]; 
ub = [50,9e16,20,60,295];
A = [];
b = [];
Aeq = [];
beq = [];

 bestx = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,[],options)
%  bestx = ga(fun,5,A,b,Aeq,beq,lb,ub,[],options)

A = bestx(1);
B = bestx(2);
C = bestx(3);
D = bestx(4);
E = bestx(5);

Pressfit = Fitt_Test(A,B,Real_Time,C,D,E);
loglog(Real_Time,Real_Pre,'*');
hold on
loglog(Real_Time,Pressfit,'r');
xlabel('Time')
ylabel('Response Data and Curve')
title('Data and Best Fitting Curve')
legend('Data','Fitted Curve')
hold off


function sse = sseval(x,Real_Time,Real_Pre)

 Qp = 10; 
Total=9e20; 
Tim= Real_Time; 
S1= 15; 
S2=50; 
Tep= 293; 
% Pre=Fitt_Test(Qp,Total,Tim,S1,S2,Tep);
Qp= x(1); 
Total=x(2); 
S1= x(3); 
S2=x(4); 
Tep= x(5);
Ton = Tim;
Pror = Fitt_Test(x(1),x(2),Ton,x(3),x(4),x(5));
 sse = sum((Real_Pre-Pror).^2);
% sse =  sqrt(mean((Pre-Pror).^2));
end



function Pressure=Fitt_Test(Qn,Tot,Time,S11,S22,Tepp)
format longg
 syms theta  P  
Q = @(v) sym(v);
p0=Q(3);
s0 = Q(1);
tao0 = Q(5)*Q(10)^(-13);
v = Q(482)*Q(10)^18;
nm=Tot;
 kb=Q(3.297623483)*Q(10)^(-24);
R = Q(19872036)*Q(10)^(-10);
T = Q(Tepp);
V = Q(436)/Q(100);
 K = Q(327)*Q(10)^17;
A = Q(3143);
Nm = A*nm;
E0=Qn  ;
 p_star = nm/(s0*v*tao0);
 Sp1=Q(S11);
 Sp2=Q(S22);
 tao = V/Sp1;
 tao1 = V/Sp2;
 
 
theta = (P/p_star)^((R*T)/E0);
d_theta = simplify(diff(theta,P));
Z = d_theta*(1/P);
F = simplify(int(Z,p0,P));

t1 = linspace(Time(1),Time(fix(length(Time)*0.21)),fix(length(Time)*0.21));
t2 = linspace(Time(fix(length(Time)*0.21)),Time(fix(length(Time)*0.5)),fix(length(Time)*0.5)-fix(length(Time)*0.21));
t3 = linspace(Time(fix(length(Time)*0.5)),Time(length(Time)),fix(length(Time))-fix(length(Time)*0.5));
syms T1
 eq1 = simplify(((T1./tao)== (log(p0/P))-((Nm/(K*V))*F)));
  eq2 = simplify(((T1./tao1)== (log(p0/P))-((Nm/(K*V))*F)));
pp1 = arrayfun(@(E)vpasolve(E,P, 1e-14), subs(eq1, T1, t1));
pp2 = arrayfun(@(E)vpasolve(E, P, 1e-14), subs(eq2, T1, t2));
pp3 = arrayfun(@(E)vpasolve(E, P, 1e-14), subs(eq2, T1, t3));
% Time=[t1';t2';t3'];
PPP=vertcat(double(real(pp1')),double(real(pp2')),double(real(pp3')));
Pressure=PPP;
end

I obtained the following outputs

[Matt J]

Well done!

[Trung]

Thank you very much for your remarks & inputs. I still need your kind help as to why the fit is still not working?

[Torsten]

@Trung

Do you think anybody can reconstruct what you do in function "Fitt_Test" to obtain the model data to be compared with "Real_pre" ?

[Trung]

Hello Mr. Torsten. I was using the this as the basis for "fitt_Test". Where You and Mr. Walter have commented on @Aiman Ahmed Q. I tried to manuplate the variables manually, "Total" seems to have the effect of moving the curves closer to each other on the Y-axis, while "Q" affects the slope.

[Trung]

I tried to the cut the data (using the first 800 data points along with "Downsample") hoping that I would be able to fut these two smaller data sections, but it didn't fit properly.

[Matt J]

There's nothing we can do to help you without understanding the prediction model.

[Trung]

Thank you for your comment. I do appreaicte all the help so far. I'm not sure if thses would help, but I have attached the main equations used in the model (Fitt_Test).

• Eq1 is equation (4), solved for Pres.
• Theta is coverage.
• P*

This is Freundlich isotherm, and the idea is to fit actual pressure vs time with the model.

[Matt J]

This looks like a 2-parameter problem and . Not sure how you got to 5 unknowns. If there are only 2 unknowns, you can probably do an exhaustive search for the optimal solution.

[Trung]

You are right. I have modified the code to optimize for only two prameters, nm (total number of molecules) and q (E, energey).

[Matt J]

Maybe you didn't understand what I meant by "do an exhaustive search". Since you now have only 2 unknown variables, you can use surf() to plot sseval as a 2D surface and see where its minimum lies. If desired, you could then run fmincon with the initial guess x0 chosen to lie where the minimum appears to be on the surf plot.

[Trung]

I will try that. Thank you very much for your kind help.