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

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.

[Learner]

Please anyone help me.

Is it possible to solve this type of problems?

Answer 1

[Star Strider]

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.

[Learner]

I am confused here. Means how can I put other 2 function outside the solver?

[Star Strider]

You don’t put the other functions outside the solver. You solve for all of them inside the ODE solver, then return only the one you want to fit, either (y(1), y(2), or y(3)), as the output of your objective function. The nonlinear parameter estimation function then does the fitting and returns the new parameters to your objective function, that integrates the ODEs each time, until the nonlinear parameter estimation function converges on a minimum-error solution.

[Learner]

Talking about the LINK you gave (lets put in my case):

You mean by following below code, I can return y(3) for fitting process, while leaving others i.e. y(1) and y(2)?

 function S = MonodKinetics1(B, t)
  function dS = DifEq(t, x)
   xdot = zeros(3,1);
   xdot(1) = bla bla;
   xdot(2) = bla bla;
   xdot(3) = bla bla;
   dS = xdot;
  end
 S = Sv(:,3);
end

[Star Strider]

Yes. Only return the one you need to fit.

[Learner]

Thanks for your help. I will try to work out, and will update you.

[Star Strider]

My pleasure!

It should work without your having to make too many changes to it.

[Learner]

I think I am doing it right!

Here what I am doing:

function S = Fit(c, t, ~)
y0 = rand(3,1); 
[T,Sv] = ode45(@DifEq, t, y0);
    function dS = DifEq(t, y)
        ydot = zeros(3,1);
        c(4) =2;
        c(5) =1.5; 
        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

Calling function is as follows:

Time = [0 1 2 3 4 5 6 7 8 9];
Sdata = [3 17 9 4 8 2 1 9 2 4];
c0 = rand(3,1) * 100;
[c,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Fit,c0,Time,Sdata);

Problem: My program is not giving any answer, but it is going busy for long long time!!

[Star Strider]

The important thing is that it’s working!

I don’t see anything in your code that could account for the slow convergence. You seem to be calling and returning everything correctly.

The problem could be your choice of ‘c0’. Do you have reason to believe that they’re so large, considering the magnitudes of c(4) and c(5)?

Interrupt lsqcurvefit, simply use:

c0 = rand(3,1);

and go again. If your regression hypersurface has a very gradual slope, starting so far away from the expected minimum could result in extremely slow convergence.

[Learner]

Yes you are right, now it is working.

I like to ask for one more favor:

1. How can I plot the experimental data points and fitted curve in one graph in my case?

2. How can I see changing of 'c values' (in workspace and in graph) as fitting process progresses? (i.e. from initial consideration to final value which program got after fitting)

[Star Strider]

My pleasure!

1. Probably the easiest way is to save the estimated parameters (a save command works well here), copy your ODE and solver system and statements to your main script workspace, then with the estimated parameters in the workspace, solve your differential equation at the values of your independent variable, xdata. Your y3 variable should now be in your workspace, so plot it against xdata. (I suggest you plot your ydata as points rather than a line.)
2. The easiest way is to print them from inside your objective function. So the first line (just above your y0 assignment) might simply to put one line as c, because that will print the values of c to the Command Window. To be a bit more formal, use fprintf(1,'\\tc(1) = %.6f\ c(2) = %.6f\ c(3) = %.6f\',c) I don’t suggest plotting during the fitting process because it slows the process significantly. If you want to write the c values serially to a mat-file, you can load the file and write code to see the fitting process later without slowing the actual fitting process itself. There may also be a display option for lsqcurvefit, so you can use the optimset function and set 'Display', 'iter-detailed' to see everything. This wouldn’t save anything, just prints to the Command Window.

[Learner]

Hi,

I did what you said on point 1.

However when I used the fitted parameter on my ODE and tried to see its fitting performance. The fitted plot of ODE is not even near to what I am fitting!!

I don't know what is the problem in that program. Can you please suggest me some solution.

Thanks in advance.

[Star Strider]

What were Rsdnrm,Rsd from lsqcurvefit? It is possible that your initial parameter estimates were sufficiently far from the actual values that the curve fit converged on a local minimum. This is a relatively common problem in nonlinear parameter estimation.

If you have an idea of what they should be, use those instead of rand that constrains them between (0,1).

If that doesn’t work (it always converges on the same parameter estimates regardless of initial estimates), check your ODE again to be sure it’s coded correctly and its initial conditions are appropriate.

[StarSign1997]

I copied and pasted the code here but an error occurs saying function value and YDATA are not equal in size. I am confused.

[Johan Sebastian Diaz Tovar]

StarSign1997 when copy the code make sure to take the transpose of Sdata, if you don't do that, you'll get that error. So: Sdata = Sdata'; Doing this you solve that.

regards.

Answer 2

[Alex Sha]

How about the results below:

Root of Mean Square Error (RMSE): 2.11836006576677

Sum of Squared Residual: 40.3870443141185

Correlation Coef. (R): 0.897567768947133

R-Square: 0.805627899852735

Adjusted R-Square: 0.740837199803646

Determination Coef. (DC): 0.805415739385938

F-Statistic: 4.2201395530756

Parameter Best Estimate

-------------------- -------------

c1 -0.607049737407678

c2 -1.63508159310052

c3 -0.00144177544567964

y1 Initial Value -8.12681195260699

y2 Initial Value -25.9407402385313