Using lsqcurvefit with ODE45 for parameter finding - I am in need of an urgent help

Hi all!

I have a set of experimental data vs. time and I want to find some parameters using lsqcurvefit function based on a series of ODEs.

My ODEs are:

Ccon(1)=k(1)*C(2)-k(2)*C(1)*C(3)-k(3)*C(1)*C(3)-k(8)*C(1);

Ccon(2)=-k(1)*C(2);

Ccon(3)=-((2*x)-z)*k(3)*C(1)*C(3)-k(9)*C(5)*C(3)+k(9)/K1*C(4)*C(7)+k(10)*C(5)*C(7)-k(10)/K2*C(6)*C(3);

Ccon(4)=x*k(3)*C(1)*C(3)+k(5)*C(1)+k(9)*C(5)*C(3)-k(9)/K1*C(4)*C(7);

Ccon(5)=x*k(2)*C(1)*C(3)+k(4)*C(1)-k(9)*C(5)*C(3)+k(9)/K1*C(4)*C(7)-k(10)*C(5)*C(7)+k(10)/K2*C(6)*C(3);

Ccon(6)=k(6)*C(1)+k(10)*C(5)*C(7)-k(10)/K2*C(6)*C(3);

Ccon(7)=(x-y+(y/2))*k(2)*C(1)*C(3)+((2*x)-z+(y/2))*k(3)*C(1)*C(3)+k(7)*C(1)+k(9)*C(5)*C(3)-k(9)/K1*C(4)*C(7)-3*k(10)*C(5)*C(7)+3*k(10)/K2*C(6)*C(3);

x,y,z,K1 and K2 values are known constants.

Ccon(i) variables are dC(i)/dt and C(i) values are the variables that change vs. time.

k(i) values are the parameters I want to find using lsqcurvefit function. I have tried to use odefit5.m function which was developed by Arthur Jutan ====> http://www.mathworks.com/matlabcentral/fileexchange/170-advstat/content/odefit5.m

This is the m file that I’ve written:

____________________________________________________________

function odefit5

%USES lsqcurvefit ode45 & XPRIME TO FIT PARAM IN ODES (STATS ASSIG 4A)

clc

clear all

% DATA points

tspan=[0 2.5 5 7.5 10 15]';

CO2data=[0 0.0297 0.0453 0.0457 0.0426 0.0461]';

yd=CO2data;

% take ave of y data

% Initial Conditions for odes

x0=[0 0.41 4.12 0 0 0 0]';

% Initial parameter guess

%disp(' enter two values for parms k1 k2')

k0=[2 .005 .000005 .2 .4 .1 .0005 .04 .006 .015]

% optimisation model is 'modode5' here--- it in turn calls ode45

% x0 is extra passed param

lb=[0 0 0 0 0 0 0 0 0 0];ub=[];% lower/upper bounds

[kfinal,resn,error,exf,out,l,jac]=lsqcurvefit(@modode5,k0,tspan,yd,lb,ub,[],x0);

kfinal

ymodel=yd+error;

figure

plot(tspan,ymodel,tspan,yd,'o')

grid

xlabel('time t');ylabel('y or x2')

title('fitted - , vs data o')

%==============

function y=modode5(k,tspan,x0)

%function y=modode5(k,tspan,x0)

%% example to fit ode parameters

% x0(init condition) is extra passed parameter

[t,x12]=ode45(@xprime,tspan,x0,[],k);

%% take 4th x12 for output y, this is the one being fitted

y=x12(:,4);

%==============

function Ccon =xprime(t,C,k)

%function xdot =xprime(t,x,flag,k)

% k is extra parameter passed from ode45

% make sure output is a col vector

K1=5.15;

K2=102.57;

x=6;

y=12;

z=6;

Ccon(1)=k(1)*C(2)-k(2)*C(1)*C(3)-k(3)*C(1)*C(3)-k(8)*C(1);

Ccon(2)=-k(1)*C(2);

Ccon(3)=-((2*x)-z)*k(3)*C(1)*C(3)-k(9)*C(5)*C(3)+k(9)/K1*C(4)*C(7)+k(10)*C(5)*C(7)-k(10)/K2*C(6)*C(3);

Ccon(4)=x*k(3)*C(1)*C(3)+k(5)*C(1)+k(9)*C(5)*C(3)-k(9)/K1*C(4)*C(7);

Ccon(5)=x*k(2)*C(1)*C(3)+k(4)*C(1)-k(9)*C(5)*C(3)+k(9)/K1*C(4)*C(7)-k(10)*C(5)*C(7)+k(10)/K2*C(6)*C(3);

Ccon(6)=k(6)*C(1)+k(10)*C(5)*C(7)-k(10)/K2*C(6)*C(3);

Ccon(7)=(x-y+(y/2))*k(2)*C(1)*C(3)+((2*x)-z+(y/2))*k(3)*C(1)*C(3)+k(7)*C(1)+k(9)*C(5)*C(3)-k(9)/K1*C(4)*C(7)-3*k(10)*C(5)*C(7)+3*k(10)/K2*C(6)*C(3);

Ccon=Ccon(:);

____________________________________________________________

I have 2 problems with MATLAB.

The 1st one is: when I tried use lower bonds, MATLAB returned an error which is “??? Error using ==> lsqncommon at 102 Levenberg-Marquardt and Gauss-Newton algorithms do not handle bound constraints and trust-region-reflective algorithm requires at least as many equations as variables; aborting.” [The interesting thing is that, I have also used lower bonds in the original m.file with adding extra k(i) values but the function has worked]

And the 2nd is: I have 5 data vs. time and the function only finds the parameters based on a single data with using this algorithm: http://www.mathworks.com/help/toolbox/optim/ug/eqn1198767322.gif

But in fact I need the optimized parameters which minimizes for all data values. I mean the objective function should be in the form of:

∑i∑t(yi,t,cal−yi,t,exp)^2 to estimate the k(i) parameters, where the index i runs through for all C(i) and the index t runs through the experimental data points.

I really need help for the problem. Do you think I should include another function file to add the second sum function for the second problem? And how can I add bonds for the k(i) values?

[KD]

Mojo, Were you able to fix the problem in your code? I am a similar problem where I have 5 sets of experimental data, with 7 odes, with 14 parameters to fit to the experimental data. I will be very grateful if you can let me have a copy of your code. I need urgent help right now. Thanks.

Answer 1

[Walter Roberson]

If you have urgent MATLAB questions, you are better off asking Mathworks Technical Support or Mathworks Consulting. This resource, MATLAB Answers, is handled by volunteers, and it can be anywhere from minutes to months before a volunteer who knows something about this topic happens upon your Question.