Asked by Maria Alvarez
on 18 Sep 2019

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

Answer by Are Mjaavatten
on 19 Sep 2019

Accepted Answer

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.

Maria Alvarez
on 20 Sep 2019

Thanks again, Are!

I made the suggested change and managed to get the simulated sea level - "simsealevel"- as curve on the plot.

However, I had to modify the "alfa" and "tau" parameters to a vector "c" to avoid some errors resulting from interactions between the lower and upper bounds of each.

I am not sure if as a result of these changes or what "tempi" is not recognised as a variable any longer.

Any pointers on where I am failing to represent this properly would be much appreciated.

Please see below for the code and errors.

CODE:

function bestalfatau = odeparam2()

%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=@(t) interp1(time,temp,t);

%Initial guess

c0= [0.2,82];

ODE_Sol= ode45(@(t,z)updateStates(t,z,c0,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];

bestc = lsqnonlin(myObjective,c0, lb, ub);

%% 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');

function f = updateStates(t,z,c,tempi)

f = (c(1).*tempi(t)-z)*(1/(c(2)));

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:

>> odeparam2

Unrecognized function or variable 'tempi'.

Error in odeparam2>objFcn (line 53)

ODE_Sol = ode45(@(t,z)updateStates(t,z,x),tempi, tSpan, z0);

Error in odeparam2>@(x)objFcn(x,time,sealevel,tSpan,z0) (line 28)

myObjective = @(x) objFcn(x,time,sealevel,tSpan,z0);

Error in lsqnonlin (line 206)

initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Error in odeparam2 (line 32)

bestc = lsqnonlin(myObjective,c0, lb, ub);

Caused by:

Failure in initial objective function evaluation. LSQNONLIN cannot continue.

Are Mjaavatten
on 21 Sep 2019

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

Maria Alvarez
on 23 Sep 2019

Thanks very much, Are!

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.