parameter identification with nlgreyest does not fit to data

Hi,

i am trying to estimate one parameter for a coupled pendelum model with nlgreyest. I am testing my implementation by generating data with a "random" parameter and trying to fit to that data by estimating the randomly selected one. Below is my code:

function Nonlinearpendulum:

function [dx,y] = NonlinearPendulum(t,x,u,m, l, g, F0, k, d0, lf, a, c, varargin)

        m1 = m;             % mass 1
        m2 = m;             % mass 2
        m3 = m;             % mass 3

        lf1  = lf;          % distance 1 from rotation axis to spring
        lf12 = lf;          % distance 2 from rotation axis to spring
        lf22 = lf;          % distance 2 from rotation axis to spring
        lf3  = lf;          % distance 3 from rotation axis to spring

        J1 = m1*l^2;        % Moment of inertia 1
        J2 = m2*l^2;        % Moment of inertia 2
        J3 = m3*l^2;        % Moment of inertia 3

        c1 = c;             % damping constant 1
        c2 = c;             % damping constant 2
        c3 = c;             % damping constant 3

        % --- Equations of motion ---

        % cartesian  postion spring
        Xf11    = lf1*sin(x(1));
        Yf11    = lf1*cos(x(1));

        Xf12    = lf12*sin(x(2))+a;
        Yf12    = lf12*cos(x(2));
        % Hinweis: XF12 und X22 identisch (ebenfalls y)
        Xf22    = Xf12;
        Yf22    = Yf12;

        Xf23    = lf3*sin(x(3))+2*a;
        Yf23    = lf3*cos(x(3));

        d11     = sqrt((Xf12-Xf11)^2+(Yf12-Yf11)^2);
        d22     = sqrt((Xf23-Xf22)^2+(Yf23-Yf22)^2);

        %
        Fm1     = -m1*g*sin(x(1));
        alp11   = atan((Xf12-Xf11)/(Yf12-Yf11))-x(1);
        Ff1     = (F0+k*(d11-d0))*sin(alp11);
        f4      = (1/J1)*(l*Fm1+lf1*Ff1-c1*x(4));

        %
        Fm2     = -m2*g*sin(x(2));
        alp12   = atan((Xf12-Xf11)/(Yf12-Yf11))-pi-x(2);
        Ff12    = (F0+k*(d11-d0))*sin(alp12);
        alp22   = atan((Xf23-Xf22)/(Yf23-Yf22))-x(2);
        Ff22    = (F0+k*(d22-d0))*sin(alp22);
        f5      = (1/J2)*(l*Fm2+lf12*Ff12+lf22*Ff22-c2*x(5));

        %
        Fm3     = -m3*g*sin(x(3));
        alp23   = atan((Xf23-Xf22)/(Yf23-Yf22))-pi-x(3);
        Ff3     = (F0+k*(d22-d0))*sin(alp23);
        f6      = (1/J3)*(l*Fm3+lf3*Ff3-c3*x(6));

        dxdt  = [x(4);x(5);x(6);f4;f5;f6];
        dx = dxdt;
        y = x(1);
    end

how i generate my data (here i picked a value of 2e-3 for c, which is the parameter i want to estimate):

clc;
clear;

m  = 0.501;     % [kg]      mass ??
l  = 13.2e-2;   % [m]       length of pendulum ??
g  = 9.81;      % [m/s²]    gravitational acceleration
F0 = 0.01;      % [N]       initial tension
k  = 7;         % [N/m]     Spring rate
d0 = 51.2e-3;   % [m]       length of spring
c  = 2e-3;    % [Nm s]    damping constant (todo, bisher schätzung)
lf = 6.8e-2;    % [m]       distance from rotation axis to spring [2.7, 6.8 ]cm
a  = 15e-2;     % [m]       distance between pendulums (= length of spring)
t = 1;
x = [deg2rad(45);deg2rad(0);0;0;0;0];
u = 2;

%[F, data] = NonlinearPendulum(t, x, u, m, l, g, F0, k, d0, lf, a, c);

expTime = 0:1e-2:130;
tic
ODE_Sol    = ode15s(@(tt,x) NonlinearPendulum(tt,x,u, m, l, g, F0, k, d0, lf, a, c),[0 1300],[deg2rad(45),deg2rad(0),0,0,0,0]);
simY       =   deval(ODE_Sol, expTime);
toc
data = simY';
data = data(:,1);

data = iddata(data,[],1e-2,'Name','Pendulum');
data.OutputName = 'Pendulum position';
data.OutputUnit = 'rad';
data.Tstart = 0;
data.TimeUnit = 's';

and here is my implementation with nlgreyest to estimate the parameter c:

m  = 0.501;     % [kg]      mass ??
l  = 13.2e-2;   % [m]       length of pendulum ??
g  = 9.81;      % [m/s²]    gravitational acceleration
F0 = 0.01;      % [N]       initial tension
k  = 7;         % [N/m]     Spring rate
d0 = 51.2e-3;   % [m]       length of spring
c  = 1e-3;    % [Nm s]    damping constant (todo, bisher schätzung)
lf = 6.8e-2;    % [m]       distance from rotation axis to spring [2.7, 6.8 ]cm
a  = 15e-2;     % [m]       distance between pendulums (= length of spring)



order         = [1 0 6];
parameters = {m, l, g, F0, k, d0, lf, a, c};
initial_states = [deg2rad(45); 0; 0; 0; 0; 0;];
Ts = 0;
nonlinear_model = idnlgrey('NonlinearPendulum',order,parameters,initial_states,Ts);
setpar(nonlinear_model,'Fixed',{true true true true true true true true false});

opt = nlgreyestOptions;
opt.Display = 'on';
opt.SearchOptions.MaxIterations = 5;
opt.SearchMethod = 'lm';
%opt.SearchOptions.Tolerance = 1e-2;

nonlinear_model = nlgreyest(data,nonlinear_model,'Display','Full',opt);

c_est = nonlinear_model.Parameters(9).Value;
[nonlinear_c_est, dnonlinear_c_est] = getpvec(nonlinear_model, 'free');

compare(data, nonlinear_model)

As i said i want to estimate c (in this case 2e-3 in the data i want to fit). In the estimation implementation i entered a different value for c, as i want to see how the estimation process is converging towards the value i used while creating the data. Unfortunately this is not the case.

It seems like the that my code and the algorithm is working, but it is not fittig the parameter to the data properly. I get really big first-order optimalities, which do not seem to decrease over the iterations. In the end i receive a c value which is very far from what it should be and the model output does also not fit at all to the data (as seen below).

I already tried different SearchMethods and tried to thighten the tolerance, but none of these worked. I never used nlgreyest before and therefore would be very glad if someone could give me a hint on how to solve my problem. Do you think this is an issue that can be solved by optimzing the nlgreyest options or is it there maybe a problem with my implementation or my understanding of the usage of nlgreyest?

Thanks a lot in advance.