How to solve multiple ODEs to fit empirical observations by optimizing multiple constants?

I have 3 ODEs and 2 parameters to be optimized to fit the ODE's to given data..

eg dA/dt = -(K1+K2)*A;

dB/dt = K1*A;

dC/dt = K2*A

where t= time and K1,K2 are constants

I have been given A,B and C vs time data..I must manipulate K1 and K2 to match the data. How do I go about doing this using optimization toolbox preferably fmincon? Please suggest a sample code if possible..

Answer 1

[Teja Muppirala]

Below is an example that does exactly what you are describing. Save it into a file and run it. First I just made some sample data, and then I fit your equations to it, getting both K1 and K2, as well as initial conditions on the data.

function fitdata
% True Values
[K1,K2,A0,B0,C0] = deal(3.5,4.2,1,2,3); 
[T,Y0] = ode45(@(t,y)[-(K1+K2)*y(1); K1*y(1); K2*y(1)],[0:0.005:1],[A0;B0;C0]);
Ymeas = Y0 + 0.1*randn(size(Y0));
figure;
plot(T,Ymeas);
hold on;
h = plot(T,nan*Ymeas,'k','linewidth',2);
minERR = Inf;
opts = optimset('fminunc');
opts.LargeScale = 'off';
Xest = fminunc(@(X)objfun(X),[0;0;0;0;0],opts);
Xest = num2cell(Xest);
[K1,K2,A0,B0,C0] = deal(Xest{:}),
legend({'A','B','C'});

          function ERR = objfun(X);
              X = num2cell(X);
              [K1,K2,A0,B0,C0] = deal(X{:});
              [T,Yest] = ode45(@(t,y)[-(K1+K2)*y(1); K1*y(1); K2*y(1)],[0:0.005:1],[A0;B0;C0]);
              ERR = sum((Ymeas(:) - Yest(:)).^2);
              if ERR < minERR
                  minERR = ERR;
                  for n = 1:3; set(h(n),'Ydata',Yest(:,n)); end
                  drawnow;
              end;
          end;
  end

[Nitin Samuel]

thank u so much!!it works so perfectly! I have to do this so for many more parameters and many more odes..can I use the same logic and generalize??thank you so much again..

Answer 2

[Ryan G]

I'm not sure how you would do this with MATLAB only but simulink design optimization would probably handle this fairly easy.

This demo shows how it can be utilized on a simple model to match a data input.

Answer 3

[Bjorn Gustavsson]

If the ODEs are that simple it should just be to integrate them analytically, then you'd simply end up with a well overdetermined least square fitting problem for K1 and K2 (perhaps you'd get A(0), B(0) and C(0) in there as unknowns too).

If the ODEs are a bit more complicated you could try a finite difference aproximation.

Answer 4

[Star Strider]

If you are looking for a way to use an ODE solver with an objective function, I have used this strategy:

function Y = objfun(B, t)     % Objective function 

      [T,Ymtx] = ode45(@DifEq, t, x0);    % Do the ODE integration 

          function dY = DifEq(t, x)       % Function ode45 integrates DifEq
               ydot(1) =  . . .;
                    . . . 
               ydot(n) =  . . .;

               dY = ydot
          end

Y = Ymtx(:,2);          % If Ymtx has more than one column, return the one you want here

end

Note that you do not have to pass the parameter vector B specifically to DifEq, since DifEq can access the B vector since it is part of objfcn.

[Nitin Samuel]

Can you please elaborate?Or can you write a sample code using same variables and parameters mentioned by me above? Im new to matlab and this will be of great help to me.Thanks

[Nitin Samuel]

Also, In the objective function I have to optimize K1 and K2 using 3 given data sets i.e A,B and C vs Time..I think your code can optimize it matching only 1 variable..