Parallel identification of ODE Parameters

Hi,

I am trying to find a clear and efficient way of performing parameter identification in a model described by non-linear ODE.

To get myself familiar with the relevant function, I am working on a toy problem that uses MATLAB own vanderpoldemo function (see code below) and the approach described by Start Strider for a similar problem here. Unfortunately, the results are pretty poor: the one parameter is estimated to be ~2 when the true value is 3. Am I doing something completely wrong?

More in general: Is there a better way to approach this problem? I would especially be interested in a way to parallelize this problem using either CPU or GPU.

Thanks,

Francesco

function funTestVDPFittingExample()

%% Create synthetic data using true parameter value u =  3;
trueMu = 3;
tspan = [0 20];
odeInitCond = [2, 0];
[data.t,data.y] = ode45(@(t,y,p) vanderpoldemo(t,y,trueMu), ...
    tspan, odeInitCond);


%% Setup the optimization
muBestGuess = 2;
lowBounds = -10;
upBounds = +10;
optFitting = optimset('Algorithm','levenberg-marquardt',...
    'TolFun',1e-8,'TolX',1e-4,'MaxIter',1e4);
muEst = lsqcurvefit(@(p,t)odeObjFun(p,t),...
    muBestGuess,data.t,data.y,lowBounds,upBounds,optFitting);

        % wrap the ODE system for compatibility with lsqcurvefit
        function yEst = odeObjFun(muTrial,t)
            ode = @(t,y) vanderpoltest(t,y);
            [t,dydt] = ode45(ode, [0 20], odeInitCond);
            % interpolation is needed to match the size of vectors across
            % different simulations (adaptive step selection might result
            % in solutions with different size as the parameter changes)
            yEst(:,1) = linterp(t,dydt(:,1),data.t);
            yEst(:,2) = linterp(t,dydt(:,2),data.t);

            % define the ODE system
            function dydt = vanderpoltest(t,y)
                dydt = [y(2);  muTrial*(1-y(1)^2)*y(2)-y(1)];
            end
         end
    
%% VALIDATE FITTING 
% generate data using estimated parameter value
dataFit = odeObjFun(muEst,data.t);      
% plot
figure(1)
plot(data.t, data.y(:,1), 'p', data.t, dataFit(:,1), '-r')
legend('data points', ['best fit \mu_{est} = ',num2str(muEst)])
end