lsqcurvefit: Function value and YDATA sizes are not equal.

Hi,

I am fitting two data from the same experiment to obtain non-linear fit parameters that describe the whole system. How can I make the value returned by the "function" the same size as YDATA used by lsqcurvefit?

function [C, D] = kinetics(k,t)                      

x0 = [1;1;1;1];
                                   
[T,Cv] = ode45(@DiffEq,t,x0);      
     function dC = DiffEq(t,x)
      
       RHOcat       = 0.23552;           
       EPS          = 0.528;            


xdot = zeros(4,1);
xdot(1) =  (RHOcat/EPS).*(-k(1).* x(1));
xdot(2) =   k(1).* x(1) - k(2).* x(2);
xdot(3) =   k(3).* x(4);
xdot(4) =   k(2) .* x(2) - k(3) .* x(4);

    dC=xdot;
 end

   C=Cv(:,1); 
   D=Cv(:,3); 
end

Y=load('Pt330.dat');
t   = Y(6:end,1);         
cNO = Y(6:end,3);         
cNO2 = Y(6:end,4);        

c = [cNO,cNO2];


tdata = Y(:,1);          
c1data = Y(:,2);         
c2data = Y(:,3);         
c3data = Y(:,4);         

k0=[1;1;1;1]; 

options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = [];
ub = [];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,c,lb,ub,options)
end

PS: Pt330.txt can be renamed as Pt330.dat to run the code.

Answer 1

[Star Strider]

The output returned by ‘kinetics’ needs to be:

C = Cv(:,[1 3]);

in this instance.

The code now works, however the fit is definitely not optimal for . I have no idea what the ‘k’ magnitudes should be, so using rand here. This gets closer ([NO] is appropriate), however there is still something wrong. (I put lower bounds on the parameters because in my experience, kinetic constants cannot be negative. If that is not true in this system, change ‘lb’ back to an empty matrix.) Also, check to be certain that ‘DifEq’ is coded correctly.

Y=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1263580/Pt330.txt');
t   = Y(6:end,1);         
cNO = Y(6:end,3);         
cNO2 = Y(6:end,4);        

c = [cNO,cNO2];


tdata = Y(:,1);          
c1data = Y(:,2);         
c2data = Y(:,3);         
c3data = Y(:,4);         

% k0=[1;1;1;1]; 
k0 = rand(4,1)*1E-2;

options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = zeros(4,1);
ub = [];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,c,lb,ub,options)

Local minimum found.

Optimization completed because the size of the gradient is less than
1e-4 times the value of the function tolerance.

k = 4×1

0.5148 0 0 0.0008

Rsdnrm = 274.3861

Rsd = 286×2

0.0100 0.0050 -0.1476 0.1029 0.1019 0.4382 0.1470 0.6181 0.1092 0.7020 0.0674 0.7532 0.0316 0.7888 0.0053 0.8188 -0.0145 0.8409 -0.0285 0.8573

ExFlg = 1

OptmInfo = struct with fields:
       iterations: 8
        funcCount: 45
         stepsize: 5.0029e+04
     cgiterations: []
    firstorderopt: 5.2194e-07
        algorithm: 'levenberg-marquardt'
          message: 'Local minimum found.↵↵Optimization completed because the size of the gradient is less than↵1e-4 times the value of the function tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 8.506067e-11,↵is less than 1e-4*options.FunctionTolerance = 1.000000e-10.'

Lmda = struct with fields:
    upper: [4×1 double]
    lower: [4×1 double]

Jmat = 572×4

0 0 0 0 -0.4185 -0.0001 0 0 -0.6285 0.0004 0 0 -0.7066 0.0015 0 0 -0.7062 0.0000 0 0 -0.6674 -0.0004 0 0 -0.6011 -0.0004 0 0 -0.5270 -0.0003 0 0 -0.4550 -0.0001 0 0 -0.3846 -0.0001 0 0

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:numel(k)
    fprintf(1, '\t\tK(%2d) = %8.5f\n', k1, k(k1))
end

		K( 1) =  0.51481
		K( 2) =  0.00000
		K( 3) =  0.00000
		K( 4) =  0.00080

tv = linspace(min(t), max(t));
Cfit = kinetics(k, tv);


figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

function C = kinetics(k,t)                      

x0 = [1;1;1;1];
                                   
[T,Cv] = ode45(@DiffEq,t,x0);      
     function dC = DiffEq(t,x)
      
       RHOcat       = 0.23552;           
       EPS          = 0.528;            


xdot = zeros(4,1);
xdot(1) =  (RHOcat/EPS).*(-k(1).* x(1));
xdot(2) =   k(1).* x(1) - k(2).* x(2);
xdot(3) =   k(3).* x(4);
xdot(4) =   k(2) .* x(2) - k(3) .* x(4);

    dC=xdot;
 end

C = Cv(:,[1 3]); 
end

I added the plot.

.

[Matterazzi]

Simply amazing! Your coding prowess is something special. Thank you once again @Star Strider :). As for the actual final results, I will revise and refine all them to match their physical meaning. Thanks a bunch!

[Star Strider]

Thank you!

As always, my pleasure!

[Matterazzi]

Even better now.

[Star Strider]

Definitely!

What are the parameter values producing that?

Did you need to change ‘DiffEq’?

[Matterazzi]

I modified the code slighty basing on my definitions of the reaction model, as:

[C, D] = kinetics(k,t)  

Such that

C = Cv(:,[1 2]); 
D = Cv(:,[3 4]); 

[Star Strider]

Thank you.

Only ‘C’ (the first output) will be ‘seen’ and used by lsqcurvefit. The second output, ‘D’, will be ignored.

This answers my question, though, since column 3 was the wrong result to return, and that was likely the reason the second output produced such a poor fit in the code I ran.

To see the entire system:

Y=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1263580/Pt330.txt');
t   = Y(6:end,1);         
cNO = Y(6:end,3);         
cNO2 = Y(6:end,4);        

c = [cNO,cNO2];


tdata = Y(:,1);          
c1data = Y(:,2);         
c2data = Y(:,3);         
c3data = Y(:,4);         

% k0=[1;1;1;1]; 
k0 = rand(4,1)*1E-2;

options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = zeros(4,1);
ub = [];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,c,lb,ub,options)

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

k = 4×1

0.5387 0.7026 0.1701 0.0049

Rsdnrm = 0.2291

Rsd = 286×2

0.0100 0.0050 -0.1576 -0.1025 0.0870 0.0487 0.1306 0.0826 0.0928 0.0500 0.0519 0.0144 0.0178 -0.0171 -0.0068 -0.0365 -0.0248 -0.0515 -0.0372 -0.0626

ExFlg = 3

OptmInfo = struct with fields:
       iterations: 14
        funcCount: 86
         stepsize: 5.0505e-04
     cgiterations: []
    firstorderopt: 8.6148e-05
        algorithm: 'levenberg-marquardt'
          message: 'Local minimum possible.↵↵lsqcurvefit stopped because the final change in the sum of squares relative to ↵its initial value is less than the value of the function tolerance.↵↵<stopping criteria details>↵↵Optimization stopped because the relative sum of squares (r) is changing↵by less than options.FunctionTolerance = 1.000000e-06.'

Lmda = struct with fields:
    upper: [4×1 double]
    lower: [4×1 double]

Jmat = 572×4

0 0 0 0 -0.4128 -0.0001 -0.0000 0 -0.6116 0.0000 0.0000 0 -0.6788 0.0000 0.0000 0 -0.6713 -0.0001 -0.0000 0 -0.6241 0.0001 0.0000 0 -0.5548 -0.0002 -0.0000 0 -0.4820 0.0002 0.0000 0 -0.4073 -0.0006 -0.0000 0 -0.3431 0.0007 0.0000 0

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:numel(k)
    fprintf(1, '\t\tK(%2d) = %8.5f\n', k1, k(k1))
end

		K( 1) =  0.53874
		K( 2) =  0.70255
		K( 3) =  0.17008
		K( 4) =  0.00490

tv = linspace(min(t), max(t));
Cfit = kinetics(k, tv);


figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

tv4 = linspace(min(t), 50);
C4 = kinetics4(k,tv4);

figure
plot(tv4, C4)
grid
xlabel('Time')
ylabel('Concentration')
legend(compose('C_%d',1:4), 'Location','best')

function C = kinetics(k,t)                      

x0 = [1;1;1;1];
                                   
[T,Cv] = ode45(@DiffEq,t,x0);      
     function dC = DiffEq(t,x)
      
       RHOcat       = 0.23552;           
       EPS          = 0.528;            


xdot = zeros(4,1);
xdot(1) =  (RHOcat/EPS).*(-k(1).* x(1));
xdot(2) =   k(1).* x(1) - k(2).* x(2);
xdot(3) =   k(3).* x(4);
xdot(4) =   k(2) .* x(2) - k(3) .* x(4);

    dC=xdot;
 end

C = Cv(:,[1 2]); 
end

function C = kinetics4(k,t)                      

x0 = [1;1;1;1];
                                   
[T,Cv] = ode45(@DiffEq,t,x0);      
     function dC = DiffEq(t,x)
      
       RHOcat       = 0.23552;           
       EPS          = 0.528;            


xdot = zeros(4,1);
xdot(1) =  (RHOcat/EPS).*(-k(1).* x(1));
xdot(2) =   k(1).* x(1) - k(2).* x(2);
xdot(3) =   k(3).* x(4);
xdot(4) =   k(2) .* x(2) - k(3) .* x(4);

    dC=xdot;
 end

C = Cv; 

end

I was just curious to see what the entire system looked like!

.

[Matterazzi]

Two more questions:

1) Curiosity killed the "bear wearing a cap" but satisfaction brought it back :)

You are ahead of me in thought. But now I wonder how I can make lsqcurvefit to 'see' the second output "D". I think it is needed for describing the system fully i.e., for NO2.

I think I am wrong in treaing lsqcurvefit as any user-defined function as function [y1,...,yN] = myfun(x1,...,xM)

2) Is there an efficient way of evaluating lsqcurvefit's goodness-of-fit wihout changing to another solver? I can surely code the usual R2 but I just checking from your vast kinetics coding experience. I read that you can combine residuals and the Jacobian to achieve this.

[Star Strider]

1) The lsqcurvefit function only needs to ‘see’ the data it is fitting. Everything else is irrelevant to it. ( I did not look closely at ‘DiffEq’ previoously. I see only three kinetic constants, so only three need to be defined in ‘x0’ and returned in ‘k’.) If you are fitting data to all four differential equations, then return all the columns, as in ‘kinetics4’. The benefit of that is that it might help define ‘k(3)’ that is not otherwise used in fitting the first two equations. If there are data for them, this simply involves re-definng ‘C’ in ‘DiffEq’.

I do not understand ‘I think I am wrong in treaing lsqcurvefit as any user-defined function ...’ so I cannot comment on it.

2) The lsqcurvefit function is the only function that I know of that will fit matrix dependent variables (there are others such as ga, however they do not return the Jacobian and other outputs necessary to calculate statistics on the fit), so it is the only option here. (I do not have the Curve Fitting Toolbox, so I do not know if it could be used for goodness-of-fit estimations involving matrix dependent variables.)

It is possible to get parameter confidence intervals using the Statistics and Machine Learning Toolbox function nlparci however not confidence intervals on the fitted regression lines, since the nlpredci function only works with nonlinear functions that return one dependent variable. The problem is that fitting a matrix of dependent variables significantly limits the options for calculating other statistics on the fit.

I have only calculated parameter statistics in kinetics fitting of matrix dependent variables, not statistics on the fit, since that is a much more complicated problem, and I am usually only interested in the parameter statistics (that tell me everything I usually need to know). I am not certain how to calculate statistics on the fit in matrix dependent variable regressions, since I have never needed to do it.