Fitting Coupled ODE's to Multiple Sets of Experimental Data

I am currently trying to fit a set of coupled kinetic differential equations to to experimentally-obtained data such that:

1. 
2. 
3. 

Where I aim to find values for , , , , and ; I already have obtained values for α, β, , and . In my MATLAB function, I have defined yy which is an n by 3 matrix containing all three sets of X vs. t, S vs. t, and P vs. t data into my function against the ODE45-determined solution. The returned solution should provide the modelled transient behavior for X, S, and P vs. t in the same order as the equations provided above.

B0 = rand(5,1)*100;
[B, Rsdnrm, Rsd, Exflg, OptmInfo, Lmda, Jmat] = lsqcurvefit(@BFERMENT,B0,tt,yy);

function Y = BFERMENT(B,t,~)
%y(1) = X, y(2) = S_lac, y(3) = P_act
%B(1) = mu_max, B(2) = K_s, B(3) = K_i, B(4) = P_m, B(5) = k_d

S0 = 41.8;
X_0 = 0.01;

X0 = [X_0, S0, 0];
tspan = [0 53];

[T,SX] = ode45(@BFERMENTFIT,tspan,X0);

    function dydt = BFERMENTFIT(t,y)

    %P production model
    alpha = 7.1328;
    beta = 0.0825;

    %S consumption model
    Yps = 0.986;
    Yxs = 0.093;

    dydt = zeros(3,1);
    dydt(1) =  ((B(1).*y(2))./((1+B(2)./y(2)).*(1+y(2)./B(3)).*(1+(y(3)./B(4)).^4))-B(5)).*y(1);
    dydt(3) =  alpha.*dydt(1)+beta.*y(1);
    dydt(2) =  (-1/Yps).*dydt(3)-(B(1)./Yxs).*y(1);   
    end

    if nargin == 2
        Y = SX(:,1);
    elseif nargin == 3
        Y = SX;
    end
end

However, whenever I try to run this script with my experimentally-obtained data, it returns an error message saying:

Error using lsqcurvefit. Function value and YDATA sizes are not equal.

Error in Data_Analysis. [B, Rsdnrm, Rsd, Exflg, OptmInfo, Lmda, Jmat] = lsqcurvefit(@BFERMENT,B0,tt,yy);

Would someone be able to help me pinpoint the issue with this? I've been racking my brain for hours trying to figure this out. It would be much appreciated.

[Star Strider]

I suspect the error may be due to:

    if nargin == 2
        Y = SX(:,1);
    elseif nargin == 3
        Y = SX;
    end

however it is not possible to determine that with the information prrovided.

If you are fitting three columns of data, the ‘BFERMENT’ function must always provide three colums of output. If you are fitting fewer columns of data, ‘BFERMENT’ must only return the matching columns.

[Nikhil Chellam]

Hi Star Strider,

Thank you for your response! I tried not having that portion in earlier because I suspected the same thing, but it still returns the same values. I can provide the lines of code I implemented earlier, though I will use some dummy numbers to make things a little bit easier.

tt = linspace(0,53,50);
yy = zeros(length(xx),3);
for i = 1:3
    yy(:,i) = rand(length(tt),1);
end

This should return two matrix quantities that have the same dimensions as what I'm trying to input. xx should be a 50x1 column vector, while yy is a 50x3 matrix. Does this provide enough information?

Thanks again!

[Star Strider]

My pleasure.

If the data you want to fit, ‘yy’, is a (50x3) matrix, and those values match what you want your objective function to fit,. return all 3 columns of the integrated diffrerential equation from your objective function.

[Nikhil Chellam]

I believe I had already tried that, doing it such that

B0 = rand(5,1)*100;
[B, Rsdnrm, Rsd, Exflg, OptmInfo, Lmda, Jmat] = lsqcurvefit(@BFERMENT,B0,tt,yy);

function Y = BFERMENT(B,t,~)
%y(1) = X, y(2) = S_lac, y(3) = P_act
%B(1) = mu_max, B(2) = K_s, B(3) = K_i, B(4) = P_m, B(5) = k_d

S0 = 41.8;
X_0 = 0.01;
X0 = [X_0, S0, 0];
tspan = [0 53];

[T,SX] = ode45(@BFERMENTFIT,tspan,X0);

    function dydt = BFERMENTFIT(t,y)
    %P production model
    alpha = 7.1328;
    beta = 0.0825;
    %S consumption model
    Yps = 0.986;
    Yxs = 0.093;
    
    dydt = zeros(3,1);
    
    dydt(1) =  ((B(1).*y(2))./((1+B(2)./y(2)).*(1+y(2)./B(3)).*(1+(y(3)./B(4)).^4))-B(5)).*y(1);
    dydt(3) =  alpha.*dydt(1)+beta.*y(1);
    dydt(2) =  (-1/Yps).*dydt(3)-(B(1)./Yxs).*y(1);   
    end
    
Y = SX;
    
end

However, it still returned with the same error. I tried just fitting it with one set of variable data, doing it such that I would pass

[B, Rsdnrm, Rsd, Exflg, OptmInfo, Lmda, Jmat] = lsqcurvefit(@BFERMENT,B0,tt,yy(:,1));

and setting

Y = SX(:,1);

At the end of the BFERMENT function to constrain the size of the output thereby ensuring the function value and YDATA sizes as equal. Nevertheless, it now appears to run ad infinitum. Is what I'm doing matching with what you suggested? Thanks again!

[Star Strider]

I have no idea what the problem could be.

[Nikhil Chellam]

After fiddling around with the code some more, it appears that the lsqcurvefit terminates when the length of the Y-output matrix does not match the size of the y-data input. As ode45 is an adaptive runge-kutta method, it looks like the number of fitted data points per iterations changes. When the length of the matrix changes, then it terminates and returns the error message, as it does not match the size of the initial matrix passed to the function. For instance, I made the function output the size of the Y-output matrix each time lsqcurvefit called BFERMENT, such that

function Y = BFERMENT(B,t)
... %everything in the previously mentioned script stays the same
    Y = SX;
disp(size(Y))
end

And after seeing the size of the Y-matrix the first time it ran, I made that the size of my test random data yy matrix to see what would happen as the script continued. The output is attached. It appears that the size of the output matrix changed, though I'm not sure how to control that. Do you have any suggestions?

[Star Strider]

That can occur if ode45 encounters a singularity and then stops its integration, so the output is not the same size as the data (dictated by the size of the ‘tspan’ vector). I have no idea what else could be causing that problem.