First order ODE parameter fit to dataset

2 次查看(过去 30 天)
Hello,
I am trying to find the optimal parameters for my first order ODE to fit a dataset.
My ODE takes the following form: ds/dt= (alfa * temp - s(t)) * (1/tau)
The parameters that I need to find are "alfa" and "tau" whereas "temp" (temperature) depends on time but I have no equation to represent it, just a dataset of its evolution over time during the same timeframe of "ds/ dt" evolution.
I would appreciate any feedback on how to get "temp" into the code as my working version is not working. Please see below the code and its errors.
Many thanks in advance for any feedback.
CODE:
%Initial data
load rcp85_expansionmidmatlab.txt;
load rcp85_temperaturemidmatlab.txt;
time= rcp85_temperaturemidmatlab(:,1)-2006+1;
temp= rcp85_temperaturemidmatlab(:,2);
sealevel=rcp85_expansionmidmatlab(:,2);
plot(time,sealevel)
%ODE information
tSpan = [1 95];
z0 = [0.0118];
tempi=temp(tSpan);
%Initial guess
alfa=0.2;%lower range from Mengel paper
tau=82;%lower range from Mengel paper
ODE_Sol= ode45(@(t,z)updateStates(t,z,alfa,tau,tempi),tSpan,z0); % Run the ODE
simsealevel = deval(ODE_Sol, time); % Evaluate the solution at the experimental time steps
hold on
plot(time, simsealevel, '-r')
%% Set up optimization
myObjective = @(x) objFcn(x,time,sealevel,tSpan,z0);
lb = [0.2,82];
ub = [0.63,1290];
bestalfatau = lsqnonlin(myObjective, alfa,tau, lb, ub);
%% Plot best result
ODE_Sol = ode45(@(t,z)updateStates(t,z,bestalfatau,tempi),tSpan, z0);
bestsealevel = deval(ODE_Sol, time);
plot(time, bestsealevel, '-g')
legend('IPCC Data','Initial Param','Best Param');
function f = updateStates(t,z,alfa,tau,tempi)
f = (alfa.*tempi-z)*(1/tau);
function cost= objFcn (x,time,sealevel,tSpan,z0)
ODE_Sol = ode45(@(t,z)updateStates(t,z,x),tempi, tSpan, z0);
simsealevel = deval(ODE_Sol, time);
cost = simsealevel-sealevel;
ERRORS:
Error using odearguments (line 95)
@(T,Z)UPDATESTATES(T,Z,ALFA,TAU,TEMPI) returns a vector of length 2, but the length of initial conditions vector is 1.
The vector returned by @(T,Z)UPDATESTATES(T,Z,ALFA,TAU,TEMPI) and the initial conditions vector must have the same number
of elements.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in odeparam1 (line 21)
ODE_Sol= ode45(@(t,z)updateStates(t,z,alfa,tau,tempi),tSpan,z0); % Run the ODE

采纳的回答

Are Mjaavatten
Are Mjaavatten 2019-9-19
The challenge is how to use your table of temperatures in your updateStates function. At a given time t, you must interpolate your table to find the correct value. Here is a simple example:
T = [22,24,18,21]; % temperature series
ts =[0,5,10,20]; % Times when T is sampled
x0 = 16;
fun = @(t,x) interp1(ts,T,t)-x;
[t,x] = ode45(fun,[0,20],x0);
plot(t,x,ts,T)
Hopefully, you will be able to use this trick to solve your ODE. The next step is to of course to find your parameters. If you have access to the optimization toolbox, try fmununc or lsqnonlin.
  5 个评论
Are Mjaavatten
Are Mjaavatten 2019-9-21
I made some modifications to make the code run with my dummy data. I have added some comments to most of my changes.
function bestc = odeparam2()
%Initial data
time = linspace(1,95,10);
temp = [20,21,23,22,27,28,23,22,23,25];
% sealevel=rcp85_expansionmidmatlab(:,2);
sealevel = sin(time);
plot(time,sealevel)
%ODE information
tSpan = [1 95];
z0 = [0.0118];
tempi=@(t) interp1(time,temp,t);
%Initial guess
c0= [0.2,82];
c0= [0.2,20];
% I find it easier to interpret the code it I define fun outside
% the call to ode45, bur this is a matter of taste
fun = @(t,z) updateStates(t,z,c0,tempi);
ODE_Sol= ode45(fun,tSpan,z0); % Run the ODE
simsealevel = deval(ODE_Sol, time); % Evaluate the solution at the experimental time steps
hold on
plot(time, simsealevel, '-r')
%% Set up optimization
myObjective = @(x) objFcn(x,time,sealevel,tempi,tSpan,z0);
% lb = [0.2,82];
% ub = [0.63,1290];
% bestc = lsqnonlin(myObjective,c0, lb, ub);
% The boundaries gave some errors that I did not bother to investigate
bestc = lsqnonlin(myObjective,c0);
%% Plot best result
ODE_Sol = ode45(@(t,z)updateStates(t,z,bestc,tempi),tSpan, z0);
bestsealevel = deval(ODE_Sol, time);
plot(time, bestsealevel, '-g')
legend('IPCC Data','Initial Param','Best Param');
end
function f = updateStates(t,z,c,tempi)
f = (c(1).*tempi(t)-z)*(1/(c(2)));
end
function cost= objFcn (x,time,sealevel,tempi,tSpan,z0)
% disp(x) % to follow the iteration process
fun = @(t,z) updateStates(t,z,x,tempi);
ODE_Sol = ode45(fun, tSpan, z0);
simsealevel = deval(ODE_Sol, time);
ds = simsealevel-sealevel;
cost = ds*ds'; % The objective must be a scalar
end
With my data, the objective was almost insenitive to tau for small values of alfa, so lsqnonln did not vary tau. With your data, it may be different. Here is how I visualised myObjective:
alfa = 0.001:0.001:0.02;
tau = 10:5:80;
[A,T] = meshgrid(alfa,tau);
r = zeros(size(A));
for i = 1:15;for j = 1:20;r(i,j) = myObjective([A(i,j),T(i,j)]);end;end
figure;surf(alfa,tau,r)
xlabel \alpha;ylabel \tau;zlabel objective

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile 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