ODE parameter optimisation to fit dataset

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!

[Pravin Jagtap]

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]

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]

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:

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.

• 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

[Maria Alvarez]

Hi Pravin,

Thank you very much for your feedback.

As I am a newbie to Matlab, I started to tackle this issue by adapting a code I found on Matlab answers that seemed similar to my case. In this code they use ODE45 to solve the ODE by setting initial values for the parameters and then LSQNONLIN to minimise the difference between the dataset and the simulated data. Please see here for the detailed explanation: https://uk.mathworks.com/matlabcentral/answers/101883-how-do-i-estimate-or-optimize-the-parameters-of-my-ode-system-in-matlab-8-1-r2013a

I might be wrong here, but if I interpret this correctly your proposed course of action assumes there is a solution for the ODE whereas I believe that might not be the case. If the latter were correct, then minimising the difference between the dataset and the simulated data would perhaps be the course of action.

[Pravin Jagtap]

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]

Thanks very much Pravin!