MATLAB Answers

Using fmincon and multistart to fit parameters of an ODE

조회 수: 2(최근 30일)
Alistair McQueen
Alistair McQueen 2021년 3월 22일
댓글: Alistair McQueen 2021년 3월 25일
Previously, I had just used fminsearch to tackle this problem.
However, especially for more complex problems (note this is rather simple) the intial guess for fminsearch is crucial, such that when its poor, convergence on a local minimum may occur.
As a result, I have been advised to use a global optimisation technique. Fminsearch is not applicable here, so I used fmincon and left the bounds 'open'.
The code here does converge to the same parameters as my fminsearch, without the optimisation but I observed a lot of integration issues and feel like the problem could be tackled better.
Was hoping to maybe get some tips/advice such that I can transition this code into my more complex/coupled ODE's.
%% Example Code - Fmincon/MultiStart
%% Section 1 - The fitting Code: Obtaining the best-fitting parameter values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%--------Fitting Set-up--------%%
function FittingCode
clc
clear all
clear workspace
p0 = [0.3,80000]; %Initial guess vector
%Syntax for error plot
%You dont need the first line - but this shows the error plot - its useful
options = optimset('Display','iter','PlotFcns',@optimplotfval);
[x,fval,exitflag,output] = fmincon(@(par)ExpData(abs(par)),p0,[],[],[],[],[0 0],[],[],options) %Initiation of fminsearch
problem=createOptimProblem('fmincon','x0',p0,'objective',@(par)ExpData(par),....
'lb',[],'ub',[]);
%Find a global solution using MultiStart with 50 iterations
ms=MultiStart;
[xmulti,errormulti]=run(ms,problem,50)
%Function where the error is calculated following ODE's (or PDE) being solved
function err=ExpData(param)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%--------Data from paper and other sources--------%%
format long
% Extracting parameters from initial guess vector
g = param(1);
K = param(2);
%Example Data - INPUT your own experimental (y-axis) data here!
cnd = [10000 18000 27000 40000 54500 70000 81000 92300 101000 104500 103900];
errcnd = [996 943 2410 2462 1363 203 1020 380 970 850 1100];
%Time data - INPUT your own time (x-axis) data here!
texp = [0 2 4 6 8 10 12 14 16 18 20];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%---------Solver set-up and data extraction---------%%
%Initial cell number - From experimental (y-axis) vector
ch0 = cnd(1);
%Solver for the ODE
[tt,xx] = ode15s(@(t,x)odeNodrug(t,x),texp,ch0);
%Data from fitting
cellC = xx(:,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%------------Error calculation------------%%
%calculate difference between simulated and experimental value
err=0; %Counter set-up
for j=1:length(texp)
err = err + ((cellC(j) - cnd(j))/cnd(j))^2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%------------Plotting------------%%
figure(1)
plot(tt,cellC, 'k', 'LineWidth',2.5)
hold on
scatter(texp, cnd, 'filled', 'k', 'Linewidth',2.5)
errorbar(texp, cnd, errcnd, 'k', 'LineStyle','none','Linewidth',2.5);
set(gca,'fontsize',32)
drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%----------Solvers----------%%
function dxdt = odeNodrug(t,x)
%The logistic growth model - NODRUG
dxdt = x(1)*g*(1 - (x(1)/K));
end
end
end

채택된 답변

Shashank Gupta
Shashank Gupta 2021년 3월 25일
Hi Alistair,
Have you tried some other global optimisation techniques such as GlobalSearch or Genetic algorithm, Check out this resource as well, it will give you an idea about how different solver works and how can it be compared.
I hope this helps or atleast gives you a headstart to explore.
Cheers.
  댓글 수: 1
Alistair McQueen
Alistair McQueen 2021년 3월 25일
Thank you - I finally managed to figure it out, using that resource funnily enough.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by