non-linear grey box System identification with Matlab

I am trying to to a non linear grey box model identification and I am using the following code. I have my measurements for the input in input vector, output vector and time stamps in time.

input = output_data(2:3,:)';
output = output_data(4:5,:)';
time = output_data(1,:)';

data = iddata(output, input, [], 'SamplingInstants', time);
data.TimeUnit = 's';

%create model
Order         = [2 2 4];               % Model orders [ny nu nx].cha
Parameters    = [1; 1; 1; 1; 1; 0.1];   % Initial parameter vector.
InitialStates = [0; 0; 0; 0];              % Initial initial states.
nlgr_m  = idnlgrey('vehicle_m', Order, Parameters, InitialStates);
setpar(nlgr_m, 'Fixed', {true true false false false false});

%Estimate the coefficients
sys = pem(data,nlgr_m, 'Display','Full', 'MaxIter', 20);

%get the parameters and the standard variation
[pvec,pvec_sd] = getpvec(sys)

I tried to use simulated input/outputs with known system parameters and the. However, the parameters that I get from this are very different from what it must be. Even when I set the initial parameter estimations It does not estimate the close parameters.

My time stamps are non-uniform which means the interval between every two sampling is not the same.

I would appreciate if anyone could help with this.

Finally, I figured out how to use nlgreyest toolbox in Matlab. Here is the code that worked for me:

M = csvread('data/all/data3.txt');
u = [M(:,5),    M(:,6)];

y = [M(:,4)* 1/10 * 3.1415/180, M(:,3) * 90/1000 * 3.1415/180 , M(:,2)];

base_elevationInit = y(1,1);
base_pitchInit = y(1,2);
base_travelInit = y(1,3);

%intial guess for the parameters
par = {-1.0000   -2.4000   -0.0943    0.1200    0.1200   -2.5000   -0.0200     0.2    2.1000   10.0000};

order = [3,2,6];  %[Ny Nu Nx]
initialStates =[base_elevationInit, base_pitchInit, base_travelInit, 0, 0, 0]';
Ts            = 0; 
m = idnlgrey('quan_model_nl',order, par, initialStates, Ts)

m.Parameters(1).Fixed = true;
m.Parameters(2).Fixed = true;
m.Parameters(8).Fixed = true;
m.Parameters(4).Fixed = true;
m.Parameters(5).Fixed = true;
m.Parameters(6).Fixed = true;
m.Parameters(9).Fixed = true;

data = iddata(y,u,0.05);

opt = nlgreyestOptions;
opt.Display = 'on';
opt.SearchOption.MaxIter = 5;

% opt.SearchMethod = 

m_est = nlgreyest(data, m, opt)

params = [m_est.Parameters(1).Value m_est.Parameters(2)

and My model function is which has to be saved in a file named quan_model_nl.m in the same folder as the previous script.

function [dx,y] = quan_model_nl(t, x, u, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, varargin)

    F = [ x(4);
           x(5);
           x(6);
           p1*cos(x(1))+ p2*sin(x(1)) + p3*x(6);
           p5*sin(x(2)) + p4*cos(x(2))+ p6*x(5);
           p7*x(6);
           ];

    G = [
                       0                   0                  ;
                       0                   0                  ;
                       0                   0                  ;                       
                       p8*cos(x(2))     p8*cos(x(2))          ;  
                       p9                -p9            ;
                       p10*sin(x(2))   p10*sin(x(2))          ;
];

    C = [
        1,0,0,0,0,0;
         0,1,0,0,0,0;
         0,0,1,0,0,0;
         ];

    dx = F + G * u';
    y = C * x ;

end