Using fmincon and multistart to fit parameters of an ODE

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

Answer 1

[Shashank Gupta]

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.

[Alistair McQueen]

Thank you - I finally managed to figure it out, using that resource funnily enough.