Fitting more than one set of data to system of ODEs

Hello! I tried to solve a kinetic model by fitting experimental data to obtain the kinetic parameters. The system consists of four time-depending variables cA, cB, cC and cD, and four ODEs containing four parameters (k1, k2, k3, k4). I want to find out the parameters k1-k4, which fit best to experimental data, consisting of several independent data sets (t1, cA1, cB1,...,t2, cA2,cB2,...).

For this purpose, I defined the ODEs:

function dc=odeFolgeParallel1(k, t, c)
dc=zeros(4,1); k=zeros(4,1);
dc(1)=-k(1)*c(1)-k(4)*c(1);
dc(2)=k(1)*c(1)-k(2)*c(2);
dc(3)=k(2)*c(2)-k(3)*c(3);
dc(4)=k(3)*c(3)+k(4)*c(1);
dc=[dc(1);dc(2);dc(3);dc(4)];

and determined the objective function using ode45

function SSE=objectiveFolgeParallel1(td, cd, k, c0)
f=@(t,c)odeFolgeParallel1(k, t, c);
[ts,cs]=ode45(f,td,c0); 
err=cd-cs; 
SSE=sum(err.^2);

After importing the experimental data

cd=[0.7 0.65 0.6;0.05 0.04 0.045;0.05 0.06 0.055;0.2 0.25 0.3]
td=[8 8.3 7.9];

and setting the initial values c0=[1;0;0;0]; k0=[0.1;0.1;0.1;0.1];

I tried to optimize my parameters ksolve using

ksolve=fminsearch(@(k)objectiveFolgeParallel1(td,cd,k,c0),k0)

getting this error

Error using odearguments (line 84)
The entries in tspan must strictly increase or decrease.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in objectiveFolgeParallel1 (line 3)
[ts,cs]=ode45(f,td,c0); %Lösung der Diff.glg.

Error in @(k)objectiveFolgeParallel1(td,cd,k,c0)

Error in fminsearch (line 189)
fv(:,1) = funfcn(x,varargin{:});

which makes sense to me because tspan has to be a timeSPAN, not specific values for the time. But how can I fix this problem, optimizing the parameters in regard of this time points? I am really stuck here and would be really grateful for some help.

Answer 1

[Torsten]

The error message is quite clear. You will have to provide td for ODE45 as td=[7.9 8 8.3].

Best wishes

Torsten.

[Teresa Schubert]

Thank you very much, I think it worked:-)

But a new problem occurs.

Error using  - 
Matrix dimensions must agree.

Error in objectiveFolgeParallel1 (line 4)
err=cd-cs; %Fehler

Error in @(k)objectiveFolgeParallel1(td,cd,k,c0)

Error in fminsearch (line 189)
fv(:,1) = funfcn(x,varargin{:});

It seems the theoretical data matrix cs and the experimental data matrix cd do nopt have the same dimension. How can I fix it so every column of cs represents the values of the solution of the ODE at the specific time of the column of cd (td)?

Answer 2

[Torsten]

err=cd'-cs;

Best wishes

Torsten.

[Teresa Schubert]

Thank you for your patience, Torsten. Fixing one error message, a new one occurs again:

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in fminsearch (line 189)
fv(:,1) = funfcn(x,varargin{:});

I am a bit stuck because I try to learn MATLAB autodidactically...

[Torsten]

SSE=sum((err()).^2);

Best wishes

Torsten.

[Teresa Schubert]

Same error again.

Is it possibly caused by the format/dimension of cd or cd somehow?

[Torsten]

Could you include the complete code you are using at the moment ?

Best wishes

Torsten.

[Teresa Schubert]

c0=[1 0 0 0]; %Defining starting contiditons
k0=[0.1 0.1 0.1 0.1];
cd=[0.650 0.0500 0.0500 0.2500; 0.6000 0.0400 0.0600 0.3000; 0.6100 0.0420 0.5800 0.2900]; %experimental data
td=[7.9;8;8.3];

Saved as odeFolgeParallel1

function dc=odeFolgeParallel1(k, t, c) %Defining ODE
  dc=zeros(4,1); k=zeros(4,1);
  dc(1)=-k(1)*c(1)-k(4)*c(1);
  dc(2)=k(1)*c(1)-k(2)*c(2);
  dc(3)=k(2)*c(2)-k(3)*c(3);
  dc(4)=k(3)*c(3)+k(4)*c(1);
  end

Saved as objectiveFolgeParallel1

function SSE=objectiveFolgeParallel1(td, cd, k, c0) %objective function
  f=@(t,c)odeFolgeParallel1(k, t, c);
  [ts,cs]=ode45(f,td,c0);  
  err=cd-cs; 
  SSE=sum((err()).^2); 
  end

Executing the Optimization

ksolve=fminsearch(@(k)objectiveFolgeParallel1(td,cd,k,c0),k0); %searching minimum of objective funtion to obtain k

Kind regards and thank you Teresa

[Torsten]

Please check the dimensions of cd and cs after the call to ode45:

[m1,n1]=size(cs);
[m2,n2]=size(cd);

Is m1 = m2 and n1 = n2 ?

Best wishes

Torsten.

[Teresa Schubert]

The problem is that I never get a result for cs because all the calculations and estimation of k have to be done to calculate cs. So it is not possible to see which dimensions cs would have.

[Torsten]

So ode45 does not calculate a solution for your starting vector k0 ?

If you try

function SSE=objectiveFolgeParallel1(td, cd, k, c0) %objective function
  f=@(t,c)odeFolgeParallel1(k, t, c);
  [ts,cs]=ode45(f,td,c0);
  [m1,n1]=size(cs);
  [m2,n2]=size(cd);
  m1
  n1
  m2 
  n2
  pause

you already get a error ?

Best wishes

Torsten.

[Teresa Schubert]

Okay, I tried it and it worked. m1=m2 (3) and n1=n2 (4), so the problem is caused by something else?

[Torsten]

Could you check the value of SSE in the same way:

function SSE=objectiveFolgeParallel1(td, cd, k, c0) %objective function
  f=@(t,c)odeFolgeParallel1(k, t, c);
  [ts,cs]=ode45(f,td,c0);
  err=cd-cs; 
  SSE=sum((err()).^2); 
  SSE
  pause

Best wishes

Torsten.

[Teresa Schubert]

Wow thank you for this hint. SSE was a row vector so it did not work with fsolve.

SSE=sum(sum(err.^2))

Calculated a skalar value for SSE and the program now runs. But for some reason, ksolve contains the same values as k0 and I do not know why.

Torsten, you are already my hero!

[Torsten]

I guess your initial conditions

c0=[1 0 0 0]; %Defining starting conditions

refer to time t=0 ?

If this is the case, you will have to set

cd=[1 0 0 0 ; 0.650 0.0500 0.0500 0.2500; 0.6000 0.0400 0.0600 0.3000; 0.6100 0.0420 0.5800 0.2900]; %experimental data

and

td=[0;7.9;8;8.3];

since you want to compare simulation and experiment at the same time instants.

Best wishes

Torsten.

[Teresa Schubert]

You are right. Unfortunately, I think there is some issue with the ODE solver as well because cs looks like

cs =

       1     0     0     0
       1     0     0     0
       1     0     0     0
       1     0     0     0

four times the initial conditions at different moments ts [0 7.9 8 8.3].

Thank you in advance!!

Teresa

[Torsten]

Of course, you will have to delete the statement

k=zeros(4,1);

in "odeFolgeParallel1" since you want to evaluate your ODEs with the actual parameters k suggested by fminsearch, not with zeros.

Best wishes

Torsten.

[Teresa Schubert]

Thank you so much!!!! It works, I can't believe it. I owe you a beer:-)

[nis]

how do we plot the curve, indictaing best fit?