MATLAB: Fitting experimental data to ODE add additional fitting parameter which should not be integrated by ODE-Solver

curve fittingexperimental dataglobal fit parametermodelode

Hey Matlab-Friends,
I want to fit some Experimental data to a First Order ODE Model and get the unknown parameters of this ODE. The fitting procedure (all fitting parameters are a part of the ode) is already working and I get my parameters as desired. Now I want to add another fitting parameter in my routine which is not part of my Model ODE. The fitting parameter is just a constant which should be added to each value of my ODE Solution, like a scaling factor.
Could anyone help me out with that? This would be very helpful. Here is the code:
 if true
% function bestc = odeParamID()
%%Create some "experimental data"
clc; clear all; clf
filename = 'Methanol_1.txt'
data_mess = dlmread(filename,'\t',0,0);
tdata = data_mess(:,1);
mdata = data_mess(:,2)/1000;
plot(tdata, mdata)
%%ODE Information
tSpan = [0 350];
z0 = 0.000000000001; % Initial values for the state
%%Initial guess
c0 = [80,84e-9,0.74,0.00001908]; % First guess at parameter value
ODE_Sol = ode15s(@(t,z)updateStates(t,z,c0), tSpan, z0); % Run the ODE
simY = deval(ODE_Sol, tdata); % Evaluate the solution at the experimental time steps
hold on
plot(tdata, simY, '-r')
%%Set up optimization
myObjective = @(x) objFcn(x, tdata, mdata,tSpan,z0);
lb = [0,84e-10,0.69,0.00001794];
ub = [180,84e-8,0.79,0.00002024];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective','OptimalityTolerance', 1.0000e-10,'FunctionTolerance',1e-10,'StepTolerance',1.0000e-10,'MaxFunEvals',10000)
bestc = lsqnonlin(myObjective, c0, lb, ub);
%%Plot best result
ODE_Sol = ode15s(@(t,z)updateStates(t,z,bestc), tSpan, z0);
bestY = deval(ODE_Sol, tdata);
plot(tdata, bestY, '-g')
legend('Exp Data','Initial Param','Best Param');
function f = updateStates(t, z, c)
g = 9.81; % m*s-2
Density = 792; %kg*m-3
SurfaceEnergy= 0.0224; %J*m-2
%c(2)r = 0.000000250; %m
%c(1)Phi= 60; %°
Viscosity= 0.000543; % kg m-2 s-1
%c(4)Surface= 0.00001;
%c(3)Porosity = 0.6;
a = 2*SurfaceEnergy/(8*Viscosity);
b = g*Density/(8*Viscosity);
q = Density;
f = a*abs(cosd(c(1)))*c(2)*q^2*c(3)^2*c(4)^2*1/z-b*c(2)^2*q*c(3)*c(4);
function cost = objFcn(x,tdata,mdata,tSpan,z0)
ODE_Sol = ode15s(@(t,z)updateStates(t,z,x), tSpan, z0);
simY = deval(ODE_Sol, tdata);
cost = simY-transpose(mdata);
  end

Best Answer

if true
% function bestc = odeParamID()
%%Create some "experimental data"
clc; clear all; clf
filename = 'Methanol_1.txt'
data_mess = dlmread(filename,'\t',0,0);
tdata = data_mess(:,1);
mdata = data_mess(:,2)/1000;
plot(tdata, mdata)
%%ODE Information
tSpan = [0 350];
z0 = 0.000000000001; % Initial values for the state
%%Initial guess
c0 = [80,84e-9,0.74,0.00001908,0]; % First guess at parameter value
ODE_Sol = ode15s(@(t,z)updateStates(t,z,c0), tSpan, z0); % Run the ODE
simY = deval(ODE_Sol, tdata); % Evaluate the solution at the experimental time steps
hold on
plot(tdata, simY, '-r')
%%Set up optimization
myObjective = @(x) objFcn(x, tdata, mdata,tSpan,z0);
lb = [0,84e-10,0.69,0.00001794,-inf];
ub = [180,84e-8,0.79,0.00002024,inf];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective','OptimalityTolerance', 1.0000e-10,'FunctionTolerance',1e-10,'StepTolerance',1.0000e-10,'MaxFunEvals',10000)
bestc = lsqnonlin(myObjective, c0, lb, ub);
%%Plot best result
ODE_Sol = ode15s(@(t,z)updateStates(t,z,bestc), tSpan, z0);
bestY = deval(ODE_Sol, tdata);
plot(tdata, bestY, '-g')
legend('Exp Data','Initial Param','Best Param');
function f = updateStates(t, z, c)
g = 9.81; % m*s-2
Density = 792; %kg*m-3
SurfaceEnergy= 0.0224; %J*m-2
%c(2)r = 0.000000250; %m
%c(1)Phi= 60; %°
Viscosity= 0.000543; % kg m-2 s-1
%c(4)Surface= 0.00001;
%c(3)Porosity = 0.6;
a = 2*SurfaceEnergy/(8*Viscosity);
b = g*Density/(8*Viscosity);
q = Density;
f = a*abs(cosd(c(1)))*c(2)*q^2*c(3)^2*c(4)^2*1/z-b*c(2)^2*q*c(3)*c(4);
function cost = objFcn(x,tdata,mdata,tSpan,z0)
ODE_Sol = ode15s(@(t,z)updateStates(t,z,x), tSpan, z0);
simY = deval(ODE_Sol, tdata);
cost = simY+x(5)-transpose(mdata);
  end
does not work ?
Best wishes
Torsten.
Related Question