ODE parameter optimisation to fit dataset

Question

Maria Alvarez 2019-9-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/480626-ode-parameter-optimisation-to-fit-dataset

评论： Maria Alvarez 2019-9-23

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 个评论
显示无隐藏无

Pravin Jagtap 2019-9-19

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 2019-9-19

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Pravin Jagtap 2019-9-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/480626-ode-parameter-optimisation-to-fit-dataset#answer_392746

编辑：Pravin Jagtap 2019-9-20

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’:

https://www.mathworks.com/help/matlab/ref/ode45.html

The possible workaround for the problem is described below:

Assumptions:

The Temperature (T) and Sea-level (S) is known as a function of Time (t).
The dS/dt can be calculated/approximated from data using any Finite difference method.
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.

This is overdetermined system and can be easily solved using ‘Least Squares method’ in MATLAB (Analytical solution exists for this kind of systems which can be coded easily). Please refer following documentation for more details: https://www.mathworks.com/help/curvefit/least-squares-fitting.html

Kind Regards

~Pravin

3 个评论
显示 1更早的评论隐藏 1更早的评论

Pravin Jagtap 2019-9-23

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 2019-9-23