Error using lsqcurvefit (line 271) Function value and YDATA sizes are not equal.

Hello everybody. I am new to MATLAB and I'm trying to find the parameters to best fit a model function having experimental points as ti and yi. The main code and functions are:

data = xlsread("Re-analyze.xlsx");

ti = data(:,6);

yi = data(:,7);

R = data(1,1);

T = data(2,1);

Me = data(3,1);

Mn0i = data(4,1);

Vs0 = data(6,1);

Mse = data(7,1);

V0 = data(8,1);

A=((36*pi)^(1/3))*((V0)^(2/3));

X0=[0.72; 28.876; 6568.86;1;1];

save variables_LpPsVb R T A Me V0 A Me Mn0i Vs0 Mse X0

function:

function X = Lp_Ps_Vb (b,t)

load variables_LpPsVb R T A V0 Me Mn0i Vs0 Mse

% % % b(1:3) = Parameters, b(4:5) = Initial Conditions

% % % MAPPING: x(1) = x, x(2) = y, b(1) = Lp, b(2) = Ps, b(3) = Vb

Ddv_Div = @(t,x,b)[((-b(1)*x(1)*A*R*T*Me*Vs0)+ b(1)*x(2)*R*T*A+b(1)*A*R*T*Mn0i*V0*Vs0-b(1)*b(3)*A*R*T*Mn0i*Vs0)/(Vs0*x(1));

(b(2)*x(1)*A*Vs0*Mse*Vs0-b(2)*x(2)*A*Vs0)/(Vs0*x(1))];

[~,X] = ode45(@(t,x) Ddv_Div(t,x,b(1:3)), t, b(4:5));

end

I get this error message:

Error using lsqcurvefit (line 271) Function value and YDATA sizes are not equal.

Unfortunately, I can not see my mistake. What should I try to solve it?

Thank you in advance!

[Mathieu NOE]

hello

the code (main section) you posted has only 13 lines and the error is about line 271

please provide the entire code otherwise it will be difficult to help you

tx

[Matt J]

@srcn The error is thrown by lsqcurvefit, but you haven't posted the code where lsqcurvefit is called - we can't see how it's being used.

[srcn]

so sorry. I wrote it in the comments below but forgot to add it here.

You can find corrected code : updated code

Answer 1

[Star Strider]

The ‘Lp_Ps_Vb’ function is going to return ‘X’ that is going to be a (numel(t) x 2) matrix, so lsqcurvefit needs to have a matching matrix of data to compare to it. (The error that was thrown can also occur when rows of data are compared to columns returned by the objective function, or the reverse, however I doubt that is the problem here.)

I am not certain what you are doing, however if you want to fit the system of differential equations to your data, the approach in Coefficient estimation for a system of coupled ODEs will likely work.

It will be necessary to adapt that code to your system. That should not be difficult.

.

[srcn]

Firstly, thank you for your answer.

I have coupled differential equations

and 3 variables which are Lp, Ps and Vb. The best fit values for these variables are Lp=0.72, Ps=28.876 and Vb= 6568.86. Other constant values (A,R,Mnoi,etc.) are obtained from the excel file I have added.

I tried to modify the code provided by you in the question : Solving Coupled Differential Equations but frankly, I was not successful.

I will try to adapt the approach you shared to my system.

[Star Strider]

What problems are you having with it?

It should be relatively straightforward to adapt.

.

[srcn]

Exactly, that should be easy to adapt but I can not see the mistakes because I don't have full knowledge of the subject. The code I wrote for the equations might be wrong or variables may not be compatible with the code. I am not sure.

That would be really great and helpful if you run the code when you have time and show the right way. I can also upload that as .mat file if you want.

Thank you very much in advance.

[Star Strider]

There are several problems with the data that it will be necessary for you to correct before I go further (if that is even possible). Also, the lsqcurvefit call is missing, so I have no idea what you are doing with it.

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/738659/Re-analyze-2.xlsx', 'VariableNamingRule','preserve')

T1 = 28×9 table
       Var1            parameters             Var3              Var4            Var5    Var6    Var7    t (min)     V (r.u.)
    __________    _____________________    __________    ___________________    ____    ____    ____    ________    ________

    {0×0 char}    {'R               ='}    8.2057e+13    {'(atm*um3/mol*K)'}    NaN     NaN     NaN            0          1 
    {0×0 char}    {'T               ='}        295.15    {'(K)'            }    NaN     NaN     NaN     0.083333    0.93461 
    {0×0 char}    {'Me           ='   }       1.8e-15    {'(mol/um3)'      }    NaN     NaN     NaN      0.16667    0.79709 
    {0×0 char}    {'Mn0i       ='     }         3e-16    {'(mol/um3)'      }    NaN     NaN     NaN         0.25    0.73734 
    {0×0 char}    {'Vbf          ='   }          0.01    {'(-)'            }    NaN     NaN     NaN      0.33333    0.72288 
    {0×0 char}    {'Vs0          ='   }    7.1001e+13    {'(um3/mol)'      }    NaN     NaN     NaN      0.41667    0.68757 
    {0×0 char}    {'Mse         ='    }       1.5e-15    {'(mol/um3)'      }    NaN     NaN     NaN          0.5    0.67377 
    {0×0 char}    {'V0            ='  }    6.5689e+05    {'(µm)'           }    NaN     NaN     NaN      0.58333    0.66694 
    {0×0 char}    {'Vend     ='       }    5.9455e+05    {'(µm)'           }    NaN     NaN     NaN      0.66667     0.6772 
    {0×0 char}    {0×0 char           }           NaN    {0×0 char         }    NaN     NaN     NaN         0.75     0.6772 
    {0×0 char}    {0×0 char           }           NaN    {0×0 char         }    NaN     NaN     NaN      0.83333    0.68757 
    {0×0 char}    {0×0 char           }           NaN    {0×0 char         }    NaN     NaN     NaN      0.91667    0.68757 
    {0×0 char}    {0×0 char           }           NaN    {0×0 char         }    NaN     NaN     NaN            1    0.69803 
    {0×0 char}    {0×0 char           }           NaN    {0×0 char         }    NaN     NaN     NaN       1.1667    0.71216 
    {0×0 char}    {0×0 char           }           NaN    {0×0 char         }    NaN     NaN     NaN       1.4167    0.71929 
    {0×0 char}    {0×0 char           }           NaN    {0×0 char         }    NaN     NaN     NaN       1.6667    0.74098 

data = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/738659/Re-analyze-2.xlsx')

data = 28×9

1.0e+13 * NaN NaN 8.2057 NaN NaN NaN NaN 0 0.0000 NaN NaN 0.0000 NaN NaN NaN NaN 0.0000 0.0000 NaN NaN 0.0000 NaN NaN NaN NaN 0.0000 0.0000 NaN NaN 0.0000 NaN NaN NaN NaN 0.0000 0.0000 NaN NaN 0.0000 NaN NaN NaN NaN 0.0000 0.0000 NaN NaN 7.1001 NaN NaN NaN NaN 0.0000 0.0000 NaN NaN 0.0000 NaN NaN NaN NaN 0.0000 0.0000 NaN NaN 0.0000 NaN NaN NaN NaN 0.0000 0.0000 NaN NaN 0.0000 NaN NaN NaN NaN 0.0000 0.0000 NaN NaN NaN NaN NaN NaN NaN 0.0000 0.0000

ti = data(:,6)

ti = 28×1

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

yi = data(:,7)

yi = 28×1

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

R = data(1,1)

R = NaN

T = data(2,1)

T = NaN

Me = data(3,1)

Me = NaN

Mn0i = data(4,1)

Mn0i = NaN

Vs0 = data(6,1)

Vs0 = NaN

Mse = data(7,1)

Mse = NaN

V0 = data(8,1)

V0 = NaN

A=((36*pi)^(1/3))*((V0)^(2/3))

A = NaN

X0=[0.72; 28.876; 6568.86;1;1]

X0 = 5×1

1.0e+03 * 0.0007 0.0289 6.5689 0.0010 0.0010

It is not possible to use xlsread with the online Run feature, however that should not be a problem with respect to reading the data.

.

[srcn]

Sorry, my bad. I updated the code.

data = [0 1

0.083333333 0.934614298

0.166666667 0.797092293

0.25 0.737335644

0.333333333 0.722878919

0.416666667 0.687565423

0.5 0.673768943

0.583333333 0.666939375

0.666666667 0.677201142

0.75 0.677201142

0.833333333 0.687565423

0.916666667 0.687565423

1 0.698034915

1.166666667 0.712161837

1.416666667 0.719293938

1.666666667 0.740980274

1.916666667 0.755676543

2.166666667 0.774319095

2.416666667 0.793265762

2.666666667 0.80093111

2.916666667 0.820307685

3.166666667 0.832081278

3.416666667 0.836031538

3.666666667 0.859993354

3.916666667 0.872142776

4.166666667 0.888520135

4.416666667 0.900938784

4.666666667 0.905101247

];

ti = data(:,1);

yi = data(:,2);

R = 8.20575E+13;

T = 295.15;

Me = 1.8E-15;

Mn0i = 3E-16;

Vs0 = 7.10015E+13;

Mse = 1.5E-15;

V0 = 656886.1948;

A=((36*pi)^(1/3))*((V0)^(2/3));

X0=[0.72; 28.876; 6568.86;1;1];

save variables R T A Me V0 A Me Mn0i Vs0 Mse ti yi

function:

function X = Lp_Ps_Vb (b,t)

load variables_LpPsVb R T A V0 Me Mn0i Vs0 Mse

% % % b(1:3) = Parameters, b(4:5) = Initial Conditions

% % % MAPPING: x(1) = x, x(2) = y, b(1) = Lp, b(2) = Ps, b(3) = Vb

Ddv_Div = @(t,x,b)[((-b(1)*x(1)*A*R*T*Me*Vs0)+ b(1)*x(2)*R*T*A+b(1)*A*R*T*Mn0i*V0*Vs0-b(1)*b(3)*A*R*T*Mn0i*Vs0)/(Vs0*x(1));

(b(2)*x(1)*A*Vs0*Mse*Vs0-b(2)*x(2)*A*Vs0)/(Vs0*x(1))];

%Ddv_Div = @(t,x,b)[-b(1)*A*R*T*(Me-x(2)/(Vs0*x(1))-Mn0i*(V0-b(3))/x(1)); Vs0*b(2)*A*(Mse-x(2)/(Vs0*x(1)))];

[~,X] = ode45(@(t,x) Ddv_Div(t,x,b(1:3)), t, b(4:5));

end

lsqcurvefit call :

B = lsqcurvefit (@Lp_Ps_Vb, X0, ti, yi, zeros(size(X0)), inf(size(X0)));

printf(1, '\n\tParameters:\n\t\tLp = %11.3E\n\t\tPs = %11.3E\n\t\tVb = %11.3E\n\tInitial Conditions:\n\t\tx = %11.3E\n\t\ty = %11.3E\n', B);

Ddv_Div = @(t,x,b)[((-b(1)*x(1)*A*R*T*Me*Vs0)+ b(1)*x(2)*R*T*A+b(1)*A*R*T*Mn0i*V0*Vs0-b(1)*b(3)*A*R*T*Mn0i*Vs0)/(Vs0*x(1));

(b(2)*x(1)*A*Vs0*Mse*Vs0-b(2)*x(2)*A*Vs0)/(Vs0*x(1))];

[Star Strider]

The problem is that the ‘Lp_Ps_Vb’ function returns two columns, however the fit is to one column of data. I created ‘Xr’ as the column to return (and arbitrarily chose the first column of ‘X’). Choose whatever column works.

This now runs, however it timed out in the online Run feature before it converged, so I was not able to estimate parameters for it here.

data = [0	1
0.083333333	0.934614298
0.166666667	0.797092293
0.25	0.737335644
0.333333333	0.722878919
0.416666667	0.687565423
0.5	0.673768943
0.583333333	0.666939375
0.666666667	0.677201142
0.75	0.677201142
0.833333333	0.687565423
0.916666667	0.687565423
1	0.698034915
1.166666667	0.712161837
1.416666667	0.719293938
1.666666667	0.740980274
1.916666667	0.755676543
2.166666667	0.774319095
2.416666667	0.793265762
2.666666667	0.80093111
2.916666667	0.820307685
3.166666667	0.832081278
3.416666667	0.836031538
3.666666667	0.859993354
3.916666667	0.872142776
4.166666667	0.888520135
4.416666667	0.900938784
4.666666667	0.905101247];

ti = data(:,1);
yi = data(:,2);
R = 8.20575E+13;
T = 295.15;
Me = 1.8E-15;
Mn0i = 3E-16;
Vs0 = 7.10015E+13;
Mse = 1.5E-15;
V0 = 656886.1948;
A=((36*pi)^(1/3))*((V0)^(2/3));
X0=[0.72; 28.876; 6568.86; 1; 1];



B = lsqcurvefit(@(b,t)Lp_Ps_Vb(b,t,R,T,Me,Mn0i,Vs0,Mse,V0,A), X0, ti, yi, zeros(size(X0)), inf(size(X0)));
fprintf(1, '\n\tParameters:\n\t\tLp = %11.3E\n\t\tPs = %11.3E\n\t\tVb = %11.3E\n\tInitial Conditions:\n\t\tx = %11.3E\n\t\ty = %11.3E\n', B);

% fprintf(1,'\tRate Constants:\n')
% for k1 = 1:length(B)
%     fprintf(1, '\t\tTheta(%d) = %14.11f\n', k1, B(k1))
% end

tv = linspace(min(ti), max(ti));
Cfit = Lp_Ps_Vb(B,tv,R,T,Me,Mn0i,Vs0,Mse,V0,A);

figure
hd = plot(ti, yi, '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, 'C_1(t)', 'Location','best')

function Xr = Lp_Ps_Vb(b,t,R,T,Me,Mn0i,Vs0,Mse,V0,A)
% load variables_LpPsVb R T A V0 Me Mn0i Vs0 Mse
      % % % b(1:3) = Parameters,  b(4:5) = Initial Conditions
      % % % MAPPING: x(1) = x,  x(2) = y,  b(1) = Lp,  b(2) = Ps,  b(3) = Vb
    Ddv_Div = @(t,x)[((-b(1)*x(1)*A*R*T*Me*Vs0)+ b(1)*x(2)*R*T*A+b(1)*A*R*T*Mn0i*V0*Vs0-b(1)*b(3)*A*R*T*Mn0i*Vs0)/(Vs0*x(1));
                                    (b(2)*x(1)*A*Vs0*Mse*Vs0-b(2)*x(2)*A*Vs0)/(Vs0*x(1))]; 
    %Ddv_Div = @(t,x,b)[-b(1)*A*R*T*(Me-x(2)/(Vs0*x(1))-Mn0i*(V0-b(3))/x(1)); Vs0*b(2)*A*(Mse-x(2)/(Vs0*x(1)))];
     
    [~,X] = ode45(@(t,x) Ddv_Div(t,x), t, b(4:5))
    Xr = X(:,1)
end

Experiment to get the correct results.

.

[srcn]

Hi Star Strider, thank you for your answer. I appreciate your help !

I've run the code but haven't been able to get the results for about an hour. I think my pc is not suitable for this analysis. I will try again on the university pc and give feedback.

Thank you very much again

[Star Strider]

It just occurred to me that this is a ‘stiff’ system (amazing what a bit of sleep can do).

data = [0	1
0.083333333	0.934614298
0.166666667	0.797092293
0.25	0.737335644
0.333333333	0.722878919
0.416666667	0.687565423
0.5	0.673768943
0.583333333	0.666939375
0.666666667	0.677201142
0.75	0.677201142
0.833333333	0.687565423
0.916666667	0.687565423
1	0.698034915
1.166666667	0.712161837
1.416666667	0.719293938
1.666666667	0.740980274
1.916666667	0.755676543
2.166666667	0.774319095
2.416666667	0.793265762
2.666666667	0.80093111
2.916666667	0.820307685
3.166666667	0.832081278
3.416666667	0.836031538
3.666666667	0.859993354
3.916666667	0.872142776
4.166666667	0.888520135
4.416666667	0.900938784
4.666666667	0.905101247];

ti = data(:,1);
yi = data(:,2);
R = 8.20575E+13;
T = 295.15;
Me = 1.8E-15;
Mn0i = 3E-16;
Vs0 = 7.10015E+13;
Mse = 1.5E-15;
V0 = 656886.1948;
A=((36*pi)^(1/3))*((V0)^(2/3));
% X0=[0.72; 28.876; 6568.86; 1; 1];
X0 = rand(5,1)*10;



B = lsqcurvefit(@(b,t)Lp_Ps_Vb(b,t,R,T,Me,Mn0i,Vs0,Mse,V0,A), X0, ti, yi, zeros(size(X0)), inf(size(X0)));

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.

fprintf(1, '\n\tParameters:\n\t\tLp = %11.6E\n\t\tPs = %11.6E\n\t\tVb = %11.6E\n\tInitial Conditions:\n\t\tx = %11.6E\n\t\ty = %11.6E\n', B);

	Parameters:
		Lp = 1.194088E-06
		Ps = 8.321523E-06
		Vb = 2.692267E+00
	Initial Conditions:
		x = 1.562840E-01
		y = 7.234029E-01

% fprintf(1,'\tRate Constants:\n')
% for k1 = 1:length(B)
%     fprintf(1, '\t\tTheta(%d) = %14.11f\n', k1, B(k1))
% end

tv = linspace(min(ti), max(ti));
Cfit = Lp_Ps_Vb(B,tv,R,T,Me,Mn0i,Vs0,Mse,V0,A);

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

function Xr = Lp_Ps_Vb(b,t,R,T,Me,Mn0i,Vs0,Mse,V0,A)
% load variables_LpPsVb R T A V0 Me Mn0i Vs0 Mse
      % % % b(1:3) = Parameters,  b(4:5) = Initial Conditions
      % % % MAPPING: x(1) = x,  x(2) = y,  b(1) = Lp,  b(2) = Ps,  b(3) = Vb
    Ddv_Div = @(t,x)[((-b(1)*x(1)*A*R*T*Me*Vs0)+ b(1)*x(2)*R*T*A+b(1)*A*R*T*Mn0i*V0*Vs0-b(1)*b(3)*A*R*T*Mn0i*Vs0)/(Vs0*x(1));
                                    (b(2)*x(1)*A*Vs0*Mse*Vs0-b(2)*x(2)*A*Vs0)/(Vs0*x(1))]; 
    %Ddv_Div = @(t,x,b)[-b(1)*A*R*T*(Me-x(2)/(Vs0*x(1))-Mn0i*(V0-b(3))/x(1)); Vs0*b(2)*A*(Mse-x(2)/(Vs0*x(1)))];
     
    [~,X] = ode15s(@(t,x) Ddv_Div(t,x), t, b(4:5));
    Xr = X(:,2);
end

I am stil not certain what column of ‘X’ should be used to fit to the data. Once that is determined, it might be worthwhile to see if the genetic algorithm (ga) function can provide a better estimate of the parameters. In my experience, ga generally works better than GlobalSearch (that I also like and is best suited for some other problems) for these sorts of problems. There are a number of others it would be worthwhile experimenting with as well.

Also, be certain thatt ‘Ddv_Div’ is coded correctly, since it does not appear to be approximating the data. (I have no idea what it is, so I cannot check it for errors.)

In any event, with the solver change, this now works and produces reasonable results in a reasonable time.

.

[srcn]

You are absolutely right. Everything is possible with a stress-free sleep (,if we can have it one day. :-) )

I would like to briefly explain what we want to obtain. We are trying to calculate water and solute permeability with the 2P-formalism which provided by Kedem-Katchalsky. These calculations are based on two main equations that I uploaded before (main equations "Ddv-Div" ).

You can find the code normally used for calculations in the attachment.(for general understanding) Normally, you will see in the code, Vb value can be calculated. However, in this case, it was desired to obtain the Vb value from two differential equations.

I will find out about the genetic algorithm you mentioned and I will carefully rewrite the differential equations.

Thank you very much once again.

[Star Strider]

I honestly don’t understand the project you’re working on. (I was an undergraduate Chemistry major, however that was almost back in the phlogiston days.) My areas of expertise very loosely (and not at all robustly) include the physics involved in the Kedem-Katchalsky equations (that I had to look up just now).

All I’m concentrating on is getting the code to work, and although tie result of the parameter estimation was not optimal (I have no control over that), I managed to do just that.

Be certain that the equations are coded correctly, and that the correct column is being returned to lsqcurvefit (or whatever optimisation function you choose) so that it has a chance of estinating the correct parameters.

That’s the best I can do.

.

[srcn]

Hi Star Strider,

It took me a long time to answer because I was doing experiments to use in this analysis.

I made a few additions and modified the code. Now it is working somehow.

I really do appreciate your wisdom on all of this. Thank you for all you have done.

[Star Strider]

As always, my pleasure!

.