idnlgrey MEX file error in the ODE file: works fine (but slow) with .m file but gets an error when working with corrensponding .c file

I want to use non linear grey-box identification on a system.

I have defined the ode equations of the system in a file called aero.m

Then I used this aero.m ode file to identify the sys. It works but it's really slow.

So I tried to write the ODE as a MEX file, by adjusting one of the mathworks examples and called this file aero_c.c

But when trying to use this .c file for the idnlgrey I get this error:

Error using iddata/nlgreyest (line 59)
The second output argument of the ODE file "aero_c" must contain 2 elements. Type
"help idnlgrey" for more information.

Here's the aero_c.c file:

/* Include libraries. */
#include "mex.h"
#include "math.h"
/* Specify the number of outputs here. */
#define NY 2

 /* State equations */
void compute_dx( double *dx, double t, double *x, double *u, double **p, const mxArray *auxvar)
{
double *dx;
/* Retrieve model parameters. */

double *M, *g, *Dm, *Dt, *Jp, *Jy, *kp, *ky, *kpy, *kyp, *dp, *dy;

M = p[0]; /* Mass of aero body */

g = p[1]; /* Gravitational acceleration */

Dm = p[2]; /* Centre of mass distance */

Dt = p[3]; /* Momentum arm */

Jp = p[4]; /* Pitch inertia */

Jy = p[5]; /* Yaw inertia */

kp = p[6]; /* Pitch thrust coefficient */

ky = p[7]; /* Yaw thrust coefficient */

kpy = p[8]; /* Pitch-yaw crossed force coeff.*/

kyp = p[9]; /* Yaw-pitch crossed force coeff.*/

dp = p[10]; /* Pitch damping coefficient*/

dy = p[11]; /* Yaw damping coefficient*/

/* x[0]: Angular pitch position */

/* x[1]: Angular pitch velocity */

/* x[2]: Angular yaw position*/

/* x[3]: Angular yaw velocity */



dx[0] = x[1];

dx[1] = (Dt[0]*kp[0]*u[0]^2+Dt[0]*kpy[0]*u[1]^2-M[0]*g[0]*Dm[0]*sin(x[0])-dp[0]*x[1]+M[0]*Dm[0]^2*x[3]^2*sin(x[0])*cos(x[0]))/Jp[0];

dx[2] = x[3];

dx[3] =  (Dt[0]*ky[0]*u[1]^2-Dt[0]*kyp[0]*u[0]^2-dy[0]*x[3])/(Jy[0]*cos(x[0]));
}

/* Output equations. */

void compute_y(double *y, double t, double *x, double **p, const mxArray *auxvar)
{

y[0] = x[0];

y[1] = x[2];

}

Here's the matlab identification code (that works fine if I put as the FileName='aero' , i.e. using the aero.m function that implement the same ODE of the aero_c.c file):

Order         = [2 2 4];             % Model orders [ny nu nx].
init          = 1e-3;
Parameters    = [M; g; Dm; Dt; Jp; Jy; init; init; init; init; init; init];   % Initial parameters values  M, g, Dm, Dt, Jp, Jy, kp, ky, kpy, kyp, dp, dy
InitialStates = x0;              % Initial initial states.
Ts            = 0;                   % Time-continuous system.

FileName      = 'aero_c'; % File describing the model structure.


% Construct the continuous-time non-linear model object
nl_m = idnlgrey(FileName, Order, Parameters, InitialStates, Ts); %initial non linear model
nl_m= setpar(nl_m, 'Fixed', {true true true true true true false false false false false false});


% Estimate the model parameters
model_non_lin = nlgreyest(Data,nl_m);

(Just for reference I also put the aero.m function with which the identification works, but too slowly for me)

function [dx, y] = aero(t, x, u, M, g, Dm, Dt, Jp, Jy, kp, ky, kpy, kyp, dp, dy, varargin)
Vp = u(1);
Vy = u(2);
x1 = x(1); %theta_pitch
x2 = x(2); %vel_pitch
x3 = x(3); %theta_yaw
x4 = x(4); %vel_yaw

% State equations
dx = [ x2; 
       (Dt*kp*Vp^2+Dt*kpy*Vy^2-M*g*Dm*sin(x1)-dp*x2+M*Dm^2*x4^2*sin(x1)*cos(x1))/Jp;
       x4;
       (Dt*ky*Vy^2-Dt*kyp*Vp^2-dy*x4/cos(x1))/(Jy*cos(x1))];
% Output equations
y = [x1;
     x3];
end

I must also say that I get the mentioned error wheter the file aero_c.c it's in the working folder or not. It's like it doesn't even exists for matlab, so I don't think it's an error directly in the aero_c.c code, but maybe I'm missing out something in the use of MEX files with this identification toolbox.

(I've also tried to run this on matlab online, regardles of the version used I always get this error)

I would be glad if someone can shed some light on this problem.

Thank you

[Walter Roberson]

Could you confirm that your .c file continues on to

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])

and so on?

[Walter Roberson]

Is there a reason why your compute_y argument list skips

    double *u,   /* Input vector (length nu). */

?

[Martin Molinaro]

For the .c file, the content it's all written in the question, I've seen that in the examples that code should work so I did not add other funtions.

I shall mention that I found an error in the ode formulation I should have use pow() for the squares, so the ode equations are ike:

dx[0] = x[1];

dx[1] = (Dt[0]*kp[0]*pow(u[0],2)+Dt[0]*kpy[0]*pow(u[1],2)-M[0]*g[0]*Dm[0]*sin(x[0])-dp[0]*x[1]+M[0]*pow(Dm[0],2)*pow(x[3],2)*sin(x[0])*cos(x[0]))/Jp[0];

dx[2] = x[3];


dx[3] =  (Dt[0]*ky[0]*pow(u[1],2)-Dt[0]*kyp[0]*pow(u[0],2)-dy[0]*x[3])/(Jy[0]*cos(x[0]));

But it's still giving me the same error.

I did not include the input *u in the compute_y function because it's not needed in the computations, I could have also not included the parameters *p since they're also not used. (it should not be a problem, at least, that's what I've understood from the mathworks examples)

To be sure I tried to put also *u in the compute_y and the rest of the function you mentioned, but it gives me the same error anyway

(to be sure the function you were referring to was this, right?:

/*DO NOT MODIFY THE CODE BELOW UNLESS YOU NEED TO PASS ADDITIONAL
   INFORMATION TO COMPUTE_DX AND COMPUTE_Y
 
   To add extra arguments to compute_dx and compute_y (e.g., size
   information), modify the definitions above and calls below.
 *-----------------------------------------------------------------------*/

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
    /* Declaration of input and output arguments. */
    double *x, *u, **p, *dx, *y, *t;
    int     i, np;
    size_t  nu, nx;
    const mxArray *auxvar = NULL; /* Cell array of additional data. */
    
    if (nrhs < 3) {
        mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidSyntax",
        "At least 3 inputs expected (t, u, x).");
    }
    
    /* Determine if auxiliary variables were passed as last input.  */
    if ((nrhs > 3) && (mxIsCell(prhs[nrhs-1]))) {
        /* Auxiliary variables were passed as input. */
        auxvar = prhs[nrhs-1];
        np = nrhs - 4; /* Number of parameters (could be 0). */
    } else {
        /* Auxiliary variables were not passed. */
        np = nrhs - 3; /* Number of parameters. */
    }
    
    /* Determine number of inputs and states. */
    nx = mxGetNumberOfElements(prhs[1]); /* Number of states. */
    nu = mxGetNumberOfElements(prhs[2]); /* Number of inputs. */
    
    /* Obtain double data pointers from mxArrays. */
    t = mxGetPr(prhs[0]);  /* Current time value (scalar). */
    x = mxGetPr(prhs[1]);  /* States at time t. */
    u = mxGetPr(prhs[2]);  /* Inputs at time t. */
    
    p = mxCalloc(np, sizeof(double*));
    for (i = 0; i < np; i++) {
        p[i] = mxGetPr(prhs[3+i]); /* Parameter arrays. */
    }
    
    /* Create matrix for the return arguments. */
    plhs[0] = mxCreateDoubleMatrix(nx, 1, mxREAL);
    plhs[1] = mxCreateDoubleMatrix(NY, 1, mxREAL);
    dx      = mxGetPr(plhs[0]); /* State derivative values. */
    y       = mxGetPr(plhs[1]); /* Output values. */
    
    /*
      Call the state and output update functions.
      
      Note: You may also pass other inputs that you might need,
      such as number of states (nx) and number of parameters (np).
      You may also omit unused inputs (such as auxvar).
      
      For example, you may want to use orders nx and nu, but not time (t)
      or auxiliary data (auxvar). You may write these functions as:
          compute_dx(dx, nx, nu, x, u, p);
          compute_y(y, nx, nu, x, u, p);
    */
    
    /* Call function for state derivative update. */
    compute_dx(dx, t[0], x, u, p, auxvar);
    
    /* Call function for output update. */
    compute_y(y, t[0], x, u, p, auxvar);
    
    /* Clean up. */
    mxFree(p);
}

Answer 1

[Martin Molinaro]

Apparently matlab had problem with seeing the .c extension, I compiled the aero_c.c file with the command:

mex aero_c.c

and this built a file: aero_c.mexw64.

I run again the code (leaving FileName='aero_c' and all like written before) and this time worked.

It did not speed up as I hoped, but it's still an improvement with respect the initial computational time.

Thank you for the suggestions!

[Walter Roberson]

MATLAB itself never automatically compiles .c files; as you noticed you have to mex the file after updating it.