MATLAB: How to fit experimental data on only one ordinary differential equations (out of multiple equations) with multiple unknown parametrs

multiple unknownsodeparameter estimation

Hi,

I am trying to estimate 3 unknown parameters that exist in 3 ordinary differential equations (ODEs) by fitting it to experimental data of only one ODE.

My equations are:

dy1/dt = (c(1)*y1.^2) + (c(2)*(y2)) – (c(3)*y3);

dy2/dt = (c(3)*y3) + (c(4)*y1);

dy3/dt = (c(1)*y1.^2) – (c(5)*y2) + (c(3)*y3);

My experimental 'x-axis data (time axis)' and 'y-axis data' for dy3/dt is as follows:

xdata = [0 1 2 3 4 5 6 7 8 9];

ydata = [3 17 9 4 8 2 1 9 2 4].

Starting guess of unknown parameters: c(1), c(2), c(3) i.e.

c0 = [1, 1, 1]; where c(4) and c(5) are known = 2 and 1.5.

Thanks in advance.

Best Answer

Definitely possible! The easiest way to describe this is to quote Monod kinetics and curve fitting.

The strategy is to put your differential equation and the ODE solver call inside the objective function you want to fit. It integrates the equations with the parameters passed to it, and returns the solved equation as the output of the objective functions. You can specify that c(4) and c(5) are known inside the objective function file, and only fit c(1)...c(3). Specify the initial conditions inside the objective function, or you can also fit them as parameters if you wish. (I’ve done all those.)

With a vector output, you can use nlinfit, lsqcurvefit, or with an extra step, fminsearch to do the fit.

Related Solutions

MATLAB: Function value and YDATA sizes are not equal

‘What am I doing wrong?’

First, not posting the actual problem, that being:

Warning: Failure at t=7.149052e+00.  Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (1.421085e-14) at time t.

so the size discrepancy is due to ode45 not being able to estimate the values for all the time elements.

Second, setting ‘c0’ to have only 3 elements when 5 are required. Fixing two of the parameters within the objective function is not appropriate. A better approach is to constrain them using the parameter bounds vectors.

Note that ‘c0’ needs to have realistic values. All the elements of ‘Sdata’ are between 1 and 17, so the last 3 elements of ‘c0’ in the slightly revised code below needs to reflect that.

I changed the code slightly to reflect this, and also to estimate the initial conditions as part of the other parameters. I leave you to experiment with it:

Time = [0 1 2 3 4 5 6 7 8 9];
Sdata = [3 17 9 4 8 2 1 9 2 4];
% Sdata = Sdata'; 
figure
plot(Time, Sdata)
grid
c0 = rand(8,1)*10;
[c,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Fit,c0,Time,Sdata);
function S = Fit(c, t, ~)
y0 = c(6:end); 
[T,Sv] = ode45(@DifEq, t, y0);
    function dS = DifEq(t, y)
        ydot = zeros(3,1);
        ydot(1) = (c(1)*y(1).^2) + (c(2)*(y(2))) - (c(3)*y(3));
        ydot(2) = (c(3)*y(3)) + (c(4)*y(1)); 
        ydot(3) = (c(1)*y(1).^2) - (c(5)*y(2)) + (c(3)*y(3));
        dS = ydot;
    end
S = Sv(:,3);
end

The problem appears to be that the differential equations may not correctly describe the system, or that the initial parameter estimates are inappropirate. Also, the system appears to be stiff, so ode15s would be an appropriate solver.

Estimating the parameters using the Global Optimization Toolbox ga (genetic algorithm) function would likely provide the best results, since lsqcurvefit requires that the initial parameter estimates are reasonably close to the actual optimal values.

MATLAB: How to fit inverse proportional function

Your function could be difficult to fit. You can use the fminsearch function to estimate the parameters:

% Define ‘x’ and ‘y’ here
% Parameter Vector: b(1) = a,  b(2) = b,  b(3) = c
yfit = @(b,x) b(1)./(x + b(2)) + b(3);  % Objective Function
CF = @(b) sum((y-yfit(b,x)).^2);        % Cost Function
b0 = rand(1,3)*10;                      % Initial Parameter Estimates
[B, fv] = fminsearch(CF, b0);           % Estimate Parameters

Note that ‘x’ and ‘y’ have to be defined in your workspace

Best Answer

Related Solutions

MATLAB: Function value and YDATA sizes are not equal

MATLAB: How to fit inverse proportional function

Related Question