Parameter estimation by fitting a system of ODEs to given data

8 Ansichten (letzte 30 Tage)
Hi, I am trying to fit a system of ODEs to my data available but I am not able to get a good fit. Can anyone help me find my mistake in the code. Thank you.
The code below is ODE function. that i am using
function S = Sfun2D(c)
% computation of an error function for an ODE model
% INPUT: b - vector of parameters
global tdata xdata x0
%c = optimvar('c',12,"LowerBound",0.1,"UpperBound",500);
%% ODE model
function dx = fan(t,x)
dx = zeros(6,1);
% c=[0.5; 0.8; 2.5; 0.9; 0.21; 0.05; 1.2; 2.3; 1; 50; 2; 0.003; 0.005; 0.05; 0.14; .00012; 1.8; 9.4; 400; 30];
% c=[0.0954; 0.074; 5.87*10^-4; 3.26*10^-6; 1.8*10^-4; 1.94*10^-4; 1.19*10-3; 3.46*10^-4; 1; 57.22; 2; 2.2; 7.27*10^-10; 0.269; 1; 1.77*10^-5; 9.4; 1.8; 1413.6; 30;];
dx(1)= (c(15)*x(5) - c(16)*(c(17)*x(2)+c(18)*x(1)+c(19)))/(c(17)*c(12)*(1+c(13)*exp(c(14)*x(1)))+c(18));
dx(2)= ((c(12)*(1+c(13)*exp(c(14)*x(1))))*(c(15)*x(5) -c(17)*(c(17)*x(2)+c(18)*x(1)+c(19))))/(c(17)*c(12)*(1+c(13)*exp(c(14)*x(1)))+c(18));
%dx(1) = ((c(15)*x(4))/(c(17)*c(12)*(1+c(13)*exp(c(14)*x(1)))+c(18)))-(eta*((c(18)+c(17)*c(12))*x(1)+((c(17)*c(12)*c(13))/c(14))*exp(c(14)*x(1))+c(17)*C+c(19))/(c(17)*c(12)*(1+c(13)*exp(c(14)*x(1)))+c(18)));
dx(3) = c(1)*x(1)+ c(20) -c(2)*x(3);
dx(4) = c(3)*(1+c(5)*x(3))-c(4)*(1+c(6)*x(3)^2)*x(4);
dx(5)= dx(1)+dx(2);
dx(6) = (c(8)*(x(3)^c(11)+c(10)^c(11)))/((c(9)*x(4)*x(3)^c(11)))- c(7)*x(6);
end
%% numerical integration set up
tspan = 0:1:max(tdata);
[tsol,xsol] = ode45(@fan,tspan,x0);
%x = linspace(0,20,50);
%y= deval(tsol,tspan)
%% plot result of the integration
figure(1)
j=0;
for i = 5:6
j=j+1;
subplot(1,2,j)
plot(tdata,xdata(:,i),'x','MarkerSize',10);
hold on
plot(tsol,xsol(:,i));
hold off
ylabel(['x(' num2str(i) ')']);
end
drawnow
%% find predicted values x(tdata)
xpred = interp1(tsol,xsol,tdata);
%% compute total error
S = 0;
for i = 1:length(tdata)
S = S + sum((xpred(i,:)-xdata(i,:)).^2);
end
end
and here is the code for the Main Run
function newparparamfit2D
% main program for fitting parameters of an ODE model to data
% the model and the error function are defined in the file Sfun2D.m
clearvars -global
global tdata xdata x0
%% data for the model
% time - value of 1st variable - value of 2nd variable
tdata(1) = 1; xdata(1,5) = 320; xdata(1,6) = 23.9;
tdata(2) = 2; xdata(2,5) = 328; xdata(2,6) = 23.8;
tdata(3) = 3; xdata(3,5) = 327.5; xdata(3,6) = 14;
tdata(4) = 4; xdata(4,5) = 319; xdata(4,6) = 15;
tdata(5) = 5; xdata(5,5) = 321; xdata(5,6) = 14;
tdata(6) = 6; xdata(6,5) = 317; xdata(6,6) = 13;
tdata(7) = 7; xdata(7,5) = 316.5; xdata(7,6) = 15;
tdata(8) = 8; xdata(8,5) = 311; xdata(8,6) = 17.5;
tdata(9) = 9; xdata(9,5) = 312; xdata(9,6) = 18;
tdata(10) = 10; xdata(10,5) = 308; xdata(10,6) = 17.9;
tdata(11) = 11; xdata(11,5) = 306; xdata(11,6) = 19;
tdata(12) = 12; xdata(12,5) = 306.3; xdata(12,6) = 20;
tdata(13) = 13; xdata(13,5) = 307; xdata(13,6) = 20;
tdata(14) = 14; xdata(14,5) = 306.3; xdata(14,6) = 20.5;
tdata(15)=15; xdata(15,5) =307; xdata(15,6) = 21.3;
tdata(16)=16; xdata(16,5) = 309; xdata(16,6) = 21.5;
tdata(17)=17; xdata(17,5) = 311; xdata(17,6) = 22;
tdata(18)=18; xdata(18,5) = 313; xdata(18,6) = 22.5;
%% initial conditions
x0= [10; 300; 2; 10.5; 320; 24];
%x0(1) = 3.5*10^5;
%x0(2) = 1;
%% initial guess of parameter values
c=[0.5; 0.8; 2.5; 0.9; 0.21; 0.05; 0.2; 2.3; 1; 3; 2; 0.003; 0.005; 0.05; 0.14; 0.00012; 1.8; 9.4; 400; 30];
% c=[0.0954; 0.074; 5.87*10^-4; 3.26*10^-6;1.8*10^-4;1.94*10^-4;1.19*10-3; 3.46*10^-4; 1;57.22;2;2.2;7.27*10^-10;0.269;1;1.77*10^-5;9.4;1.8;1413.6;30];
%b(1) = 0.01;
%b(2) = 0.2;
%% minimization step
[cmin, Smin] = fminsearch(@Sfun2D,c);
disp('Estimated parameters c(i):');
disp(cmin)
disp('Smallest value of the error S:');
disp(Smin)
end
  3 Kommentare
deva siva sai murari Kanumoori
Hello Star Strider!! thank you very much for replying. Firstly, thanks a lot for all your contributions. I think its only becasue of your answers i was able to write thins code. I had followed the methods written by you from previously asked questions and written the following code. The fminsearch works well but its not able to give me proper fit. I am not able to figure out why.
Star Strider
Star Strider am 22 Okt. 2021
My pleasure!
First — PLEASE DO NOT USE GLOBAL VARIABLES! This is especially true for optimisation problems.
Second — There are a large number of other problems.
  • The time vector if the original data should be an argument to ‘Sfun2D’ so that the simulation times match exactly with the data times. Use that as the ‘tspan’ input to the numeric differential equation integrator.
  • I always let the optimisation routine estimate the initial conditions of the differential equations as parameters, and append them as the last elements of the parameter vector. Re-define the parameter vector inside the objective function to accommodate that change.
  • The fminsearch function is not well-suited to optimise more than a few parameters (if I remember correctly, 5 is esentially its upper limit), so estimating 22 of them (my visual count) is clearly beyond its limits, and will be difficult for any optimisation routine to fit. (I would suggest using the ga function for this problem, then be prepared to wait several hours for it to converge. Set the number of generations at 10000 to give it time to explore.)
  • Rather than creating ‘xdata’ with only columns 5 and 6 filled, just fit the last 2 columns of the differential equaiton to the 2 column data matrix.
The code is difficult for me to follow, and therefore impossible for me to debug.
.

Melden Sie sich an, um zu kommentieren.

Antworten (0)

Produkte


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by