Optimizing ODE parameter to make solution fit empirical observations

Hello everyone,

I am trying to estimate the k parameter of a simple linear model (dydt = k) based on experimental observations using fmincon. The fmincon function minimizes my objective function that describes the error of the fit of the ODE solution for a certain k value and the observed values. The objective function in turn cals the ode45 solver for solving the model.

Note that I know I could also fit the solution of the ODE to the data or use polyfit to obtain parameter estimates. However, I want to take the apprach of solving the ODE using ode45 to gain experience for future applications.

The code seems to run, but however I am not obtaining a good fit. Does someone know what the problem might be or could someone give me advise to solve it? Please find the script below:

function parameters = myOptimization(measuredValues,Time)

measuredValues = [2; 4; 6; 10; 11; 15; 17; 20; 23; 25];     % Ficitve observations
Time= [1:1:10];                                             % Time corresponding with the observations 

hold on;
plot(Time,measuredValues);
h = plot(Time,nan*measuredValues,'r');
set(h,'tag','solution');

initialConditions = [3 2];
lb = [-10 -10];
ub = [10 10];
F = @(initials) COST(initials,Time,measuredValues);
parameters = fmincon(F,initialConditions,[],[],[],[],lb,ub); %  fmincon used as optimizer
legend({'Measured','Fitted'});
disp(['fmincon: parameters = ' num2str(parameters)]);

function COST = COST(initialConditions,Time,measuredValues)

y0 = initialConditions(1);
k0 = initialConditions(2);

% The cost function that calls the ODE solver.
[t,y] = ode15s(@myModel,Time,y0);

delta= (y - measuredValues).^2;
COST = delta'*delta;

%COST = sum((P - measuredValues).^2)

h = findobj('tag','solution');
set(h,'ydata',y);
title(['y0 = ' num2str(y0) '   ' 'k = ' num2str(k0)]);
drawnow;

function dydt = myModel(t,k)
dydt = k;

Answer 1

[Torsten]

[t,y] = ode15s(@(t,y)myModel(t,y,k0),Time,y0);

function dydt = myModel(t,y,k)
dydt = k;

Best wishes

Torsten.

[Niels Bessemans]

Thank you very much for your quick and useful reply!

Best wishes, Niels Bessemans

Answer 2

[Niels Bessemans]

Hi everyone,

As continuation of my previous exercise posted above, I am now trying to fit an exponential decay model to experimental data. The model is given by:

function dPdt = myModel(t,P,dP,k)
V = 400;    % m3
T = 284.15; % K
R = 8.314;  % J/molK
dPdt = ((-k*R*T)/V)*dP;

Ans my solver call looks as follows:

[t,P] = ode15s(@(t,P)myModel(t,P,dP,k0),Time,P0);

Initial conditions remained the same as in the example above. When I run the code I get the error :

">> myLeak
Error using odearguments (line 95)
@(T,P)MYMODEL(T,P,DP,K0) returns a vector of length 15, but the length of initial conditions vector is 1. The vector returned by
@(T,P)MYMODEL(T,P,DP,K0) and the initial conditions vector must have the same number of elements.

Does someone have an idea about the cause or how I could solve this problem?

[Torsten]

dP and k must be scalars in the evaluation of dPdt. But either dP or k is a vector of length 15.

Best wishes

Torsten.

Answer 3

[Niels Bessemans]

Dear Torsten,

That is right, dP is a vector of length 15 representing the measured pressure differences between measurement object and environment, which is used as input for my model to be fit.

How could I "tell" Matlab to use one element of that vector for each experimental observation?

Kind regards, Niels

[Torsten]

[t,P] = ode15s(@(t,P)myModel(t,P,Time,dP,k0),Time,P0);

function dPdt = myModel(t,P,t_array,dP_array,k)
dP = interp1(t_array,dP_array,t);
V = 400;    % m3
T = 284.15; % K
R = 8.314;  % J/molK
dPdt = ((-k*R*T)/V)*dP;

And please don't open new "Answer" threads if you want to place a "Comment".

Best wishes

Torsten.

[Niels Bessemans]

Dear Torsten,

Thanks! I'll keep it in mind to post as comments.

Kind regards, Niels