Parameter estimation by fitting a system of ODEs to given data
10 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
Star Strider
2021 年 10 月 22 日
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.
.
回答 (0 件)
参考
カテゴリ
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!