Optimizing ODE parameter to make solution fit empirical observations

2 次查看（过去 30 天）
Niels Bessemans
Niels Bessemans 2018-5-4
评论： Niels Bessemans 2018-5-9
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;

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2018-5-4
[t,y] = ode15s(@(t,y)myModel(t,y,k0),Time,y0);
function dydt = myModel(t,y,k)
dydt = k;
Best wishes
Torsten.

更多回答（2 个）

Niels Bessemans
Niels Bessemans 2018-5-9
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?
  1 个评论
Torsten
Torsten 2018-5-9
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.

请先登录，再进行评论。

Niels Bessemans
Niels Bessemans 2018-5-9
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
  2 个评论
Torsten
Torsten 2018-5-9
[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
Niels Bessemans 2018-5-9
Dear Torsten,
Thanks! I'll keep it in mind to post as comments.
Kind regards, Niels

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by