Fitting global parameters using fminsearch with nested ode

Hello,

I have a problem fitting a parameter to experimental data. The model equations build a system of ODE´s. The simplified version looks like this (A and B are model parameters, C is a constant):

dx1_dt = -A*(x1-x2)
dx2_dt = A*(x1-x2) - B*(x2-C)

The task is to estimate the optimal value of A and B for both, x1(t) and x2(t).

Currently I use fminsearch with a nested ode45 in a loop. I calculate the RMSE for x1 and x2, then add them together to a global RMSE and give it back to fminsearch as quality criteria. It works quite well if you´re not in a hurry.

But I believe theres a better way to do it.

Please consider that I don´t have access to any toolboxes. I am aware of fminsearch´s drawbacks. More generally I am wondering how the algorithm should look like.

Thanks in advance.

mulm

[Torsten]

There is no easier way to proceed as the one you describe above.

You can speed up the calculations by solving your system of ODEs analytically (which is at least possible for your simplified version from above) and/or to use the optimization toolbox. But I think both options are not possible for you.

For a guideline, take a look at

http://www.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

Best wishes

Torsten.

[Matthias Müller]

Thanks for your response, Torsten. At least my initial guess was adequate..

Regards,

mulm

Answer 1

[Star Strider]

I refer you to: Monod kinetics and curve fitting. I used lsqcurvefit here, but you can use fminsearch with slight modification of the code.

[Matthias Müller]

It does not look too trivial to implement this algorithm using fminsearch but I will try. I was hoping for a "simple" solution ;)

Thanks anyway.

[Star Strider]

My pleasure!

You can create a least-squares estimator using fminsearch fairly easily:

x = ...;                                            % Independent Variable Data                     
y = ...;                                            % Dependent Variable Data
obj_fcn = @(B,xdata) ...;                           % Objective Function (Incorporating Your ODE)

SSECF = @(B) sum((y - obj_fcn(B,xdata)).^2);        % Sum-Squared-Error Cost Function
B_init = [ ... ];                                   % Initial Parameter Estimate Vector
[B_hat,SSE] = fminsearch(SSECF, B_init);            % Estimate Parameters (‘SSE’ Is The Sum-Squared-Error At Convergence)

y_fit = obj_fcn(B_hat, x);

figure(1)
plot(x, y, 'bp')                                    % Plot Data
hold on
plot(x, y_fit, '-r')                                % Plot Fitted Data
hold of
grid

I just sketched it out here, so you will likely have to modify it a bit to make it work with your code.