Borrar filtros
Borrar filtros

ODE parameter optimisation to fit dataset

9 visualizaciones (últimos 30 días)
Maria Alvarez
Maria Alvarez el 16 de Sept. de 2019
Comentada: Maria Alvarez el 23 de Sept. de 2019
Hello,
I am trying to find the optimal values for two parameters of a first order differential equation but am having some errors that are difficult to bypass. I would be very grateful if anyone could point to sources or give me feedback on this. Please see below for the code and errors:
function bestalfatau = odeparam1()
%Initial data
load rcp85_expansionmidmatlab.txt;
load rcp85_temperaturemidmatlab.txt;
time= rcp85_temperaturemidmatlab(:,1);
temp= rcp85_temperaturemidmatlab(:,2);
sealevel=rcp85_expansionmidmatlab(:,2);
plot(time,sealevel)
%ODE information
tSpan = [2006:1:2100];
z0 = 0.0118;
%Initial guess
alfa=0.2;
tau=82;
ODE_Sol= ode45(@(t,z)updateStates(t,z,alfa,tau), 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), 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)
f = (alfa.*temp-z)*(1/tau);
function cost= objFcn (x,time,sealevel,tSpan,z0)
ODE_Sol = ode45(@(t,z)updateStates(t,z,x), tSpan, z0);
simsealevel = deval(ODE_Sol, time);
cost = simsealevel-sealevel;
ERRORS:
>> odeparam1
Unrecognized function or variable 'temp'.
Error in odeparam1>updateStates (line 49)
f = (alfa.*temp-z)*(1/tau);
Error in odeparam1>@(t,z)updateStates(t,z,alfa,tau) (line 19)
ODE_Sol= ode45(@(t,z)updateStates(t,z,alfa,tau), tSpan, z0); % Run the ODE
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in odeparam1 (line 19)
ODE_Sol= ode45(@(t,z)updateStates(t,z,alfa,tau), tSpan, z0); % Run the ODE
Many thanks in advance!
  2 comentarios
Pravin Jagtap
Pravin Jagtap el 19 de Sept. de 2019
Hello Maria,
It will be helpful if you provide the dummy/sample dataset and input parameters which you used in your code (text files) so that same errors can be reproduced. Most of the errors are self explanatory.
Maria Alvarez
Maria Alvarez el 19 de Sept. de 2019
Thank you Pravin. To give you a little bit of background, 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" (by minimising the difference between the ode calculation and the dataset) 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.
My datasets (please find attached) consist of an evolution of "temp" over time (2006- 2100) as well as of variable "s" over time (2006- 2100).
In the code I have made some changes to include "temp" as a variable. Unlike "tau" and "alfa", I think I need to include the different values of "temp" as the calculation moves through "tSpan". 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!
function bestalfatau = odeparam1()
%Initial data
load rcp85_expansionmidmatlab.txt;
load rcp85_temperaturemidmatlab.txt;
time= rcp85_temperaturemidmatlab(:,1)-2006+1;%modifying the vector from [2006 2100] to [1 95] to match tSpan
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;
tau=82;
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

Iniciar sesión para comentar.

Respuesta aceptada

Pravin Jagtap
Pravin Jagtap el 20 de Sept. de 2019
Editada: Pravin Jagtap el 20 de Sept. de 2019
Hello Maria,
The procedure followed in the code is not consistent with the objectives mentioned in the description of a problem. The ‘ode45’ is used to solve the ‘Ordinary differential equation’ not to get the parameters. Please refer following documentation for more details on ‘ode45’:
The possible workaround for the problem is described below:
  • Assumptions:
  1. The Temperature (T) and Sea-level (S) is known as a function of Time (t).
  2. The dS/dt can be calculated/approximated from data using any Finite difference method.
  3. The given ODE governs the physical phenomena described in the problem.
  • Remodify the ODE as bellow:
dS/dt = (alpha * T - S) * (1/tau)
S = alpha*T - tau*dS/dt
In the above equation, S, T and dS/dt is known and System can be expressed as linear combination of ‘alpha’ and ‘tau’. Therefor, this can be transformed as ‘Ax=b’ and solved for x, where x is [alpha tau], A = [ T dS/dt ] and b = S.
Kind Regards
~Pravin
  3 comentarios
Pravin Jagtap
Pravin Jagtap el 23 de Sept. de 2019
Hello Maria,
If I understand correctly, the MATLAB answer case which you shared has well defined objective function and optimization needs the solution of system of ODEs. Also, the concern raised about nonexistence of solution for ODE is not valid because 'ode45' approximate the solution numerically. Moreover, in the approach of "minimizing the difference between the dataset and the simulated data" assumes the approximate solution of ODE (simulation results). The solution proposed earlier is just a simple possible workaround. I would recommend trying both methods and check which one works. (I suspect that both may converge to same solution because in both the methods follow error minimization). Please refer following error free code:
function [bestalfa, besttau] = working()
%% Initial data
load One.txt;
load Two.txt;
time = Two(:,1)-2006+1; %modifying the vector from [2006 2100] to [1 95] to match tSpan
temp = Two(:,2);
sealevel = One(:,2);
plot(time,sealevel,'dm')
%% ODE information
tSpan = [1 95];
z0 = 0.0118;
%% Initial guess
alfa1= 0.2;
tau1 = 82;
ODE_Sol= ode45(@(t,z)updateStates(z,alfa1,tau1,temp,t),tSpan,z0); % Run the ODE
simsealevel = deval(ODE_Sol, time); % Evaluate the solution at the experimental time steps
hold on
plot(time, simsealevel, '-b')
%% Set up optimization
myObjective = @(x) objFcn(x,time,sealevel,tSpan,z0);
lb = [0.2,82];
ub = [0.63,1290];
best = lsqnonlin(myObjective, [alfa1 tau1], lb, ub)
bestalfa = best(1)
besttau = best(2)
%% Plot best result
ODE_Sol = ode45(@(t,z)updateStates(z, bestalfa, besttau,temp,t),tSpan, z0);
bestsealevel = deval(ODE_Sol, time);
plot(time, bestsealevel, '*g')
legend('IPCC Data','Initial Param','Best Param');
%%
function f = updateStates(z,alfa,tau,temp,t)
%t
T = interp1(1:95,temp,t);
f = (alfa*T-z)*(1/tau);
end
%%
function cost = objFcn (x,time,sealevel,tSpan,z0)
ODE_Sol = ode45(@(t,z)updateStates(z,alfa1,tau1,temp,t),tSpan,z0);
simsealevel = deval(ODE_Sol, time);
cost = simsealevel-sealevel;
end
end
I tried with different initialization value for ‘alpha’ and ‘tau’ but not happy with the solution. Please make sure your objective is consistent with code.
Maria Alvarez
Maria Alvarez el 23 de Sept. de 2019
Thanks very much Pravin!

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by