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!
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
采纳的回答
Pravin Jagtap
2019-9-20
编辑：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’:
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 个评论
Maria Alvarez
2019-9-20
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Systems of Nonlinear Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
