Getting segmentation fault when trying to run an Ipopt optimization with linear Mumps

[Alfoesp97]

Hi, i'm pretty new to mathematical optimization and i don't have much confidence in C so i'm wasting so much time to figure out how to run an optimization with Ipopt in C. First of all i've compiled and installed all the libraries and i'm using the functions of IpStdCInterface.h. When i try to run my code, I'm getting a segmentation fault on the line "IpoptSolve(nlp, x, NULL, &obj, mult_g, mult_x_L, mult_x_U, &user_data);", after some hours of debugging i found out that the error line that gives me the segmentation fault is in IpTNLPAdapter.cpp and it is the following line of code:"if( c_col != -1 && c_row != -1 )" but i can't figure out how to solve this. I tried different options regarding the jacobians computation but i'm getting the same error over and over. Below is the main of my code:

MyUserData user_data;
    user_data.ds=(double*)malloc((size_ds)*sizeof(double));
    user_data.a=(double*)malloc((size_a)*sizeof(double));
    user_data.b=(double*)malloc((size_b)*sizeof(double));
    user_data.mrif=mrif;
    user_data.m_map=m_map;
    user_data.Kref=Kref;
    user_data.tmax=tmax;
    
    for(int i=0;i<size_ds-1;i++){
        user_data.ds[i]=array_ds[i];
    }
    for (int i=0;i<size_a-1;i++){
        user_data.a[i]=a[i];
    }
    for (int i=0;i<size_b-1;i++){
        user_data.b[i]=b[i];
    }
    for (int k=0;k<3;k++){
        for(int j=0;j<3;j++){
            user_data.fmapcoeff[k][j]=fmapcoeff[k][j];
        }
    }
    free(json_data);
    cJSON_Delete(json);
    ipindex n = size_guess+1;                      /* number of variables */
    ipindex m = size_guess+2;                      /* number of constraints */
    ipindex nele_jac=0;                    /* number of nonzeros in the Jacobian of the constraints */
    ipindex nele_hess=0;                   /* number of nonzeros in the Hessian of the Lagrangian (lower or upper triangular part only) */
    ipindex index_style=0;                 /* indexing style for matrices */
    ipnumber* x_L = (ipnumber*)malloc(sizeof(ipnumber) * n);                /* lower bounds on x */
    ipnumber* x_U = (ipnumber*)malloc(sizeof(ipnumber) * n);
    
    for(int i=0;i<size_ub;i++){
        x_U[i]=ub[i];
        x_L[i]=lb[i];
    }
    /*for(int i=0;i<size_ub;i++){
        printf("%f\n ", x_U[i]);
        printf("%f\n ", x_L[i]);
    }  */              /* upper bounds on x */
    ipnumber* g_L =(ipnumber*)malloc(sizeof(ipnumber) * m);                   /* lower bounds on g */
    ipnumber* g_U = (ipnumber*)malloc(sizeof(ipnumber) * m);   
    for(int i=0;i<n;i++){
        g_L[i]=-1e-20;
        g_U[i]=0.0;
    }
    g_L[n]=user_data.tmax;
    g_U[n]=user_data.tmax;               /* upper bounds on g */
    IpoptProblem nlp = NULL;             /* IpoptProblem */
    enum ApplicationReturnStatus status; /* Solve return code */
    ipnumber* x = (ipnumber*)malloc(sizeof(ipnumber) * n);
    for(int i=0;i<n;i++){
        if(i==0)
            x[i]=k0;
        else
            x[i]=guess[i-1];
    }
                   /* starting point and solution vector */
    ipnumber* mult_g = NULL;             /* constraint multipliers at the solution */
    ipnumber* mult_x_L = NULL;           /* lower bound multipliers at the solution */
    ipnumber* mult_x_U = NULL; 
    mult_g = (ipnumber*) malloc(sizeof(ipnumber) * m);
    mult_x_L = (ipnumber*) malloc(sizeof(ipnumber) * n);
    mult_x_U = (ipnumber*) malloc(sizeof(ipnumber) * n);          /* upper bound multipliers at the solution */
    ipnumber obj;                        /* objective value */
             /* our user data for the function evaluations */
    ipindex i;                           /* generic counter */
    // Function to calculate the difference between two elements in s_sweep
    
    
    
    nlp = CreateIpoptProblem(n, x_L, x_U, m, g_L, g_U, 1, 0, index_style,
                            &eval_f, &eval_g,&eval_grad_f, &eval_jac_g,&eval_h);
    if (!x_L || !x_U || !g_L || !g_U) {
    fprintf(stderr, "Memory allocation failed.\n");
    return -1;
    }
    if (nlp == NULL) {
    fprintf(stderr, "Failed to create Ipopt problem.\n");
    return -1;
    }
    

    free(x_L);
    free(x_U);
    free(g_L);
    free(g_U);
    
    //(nlp, "mu_strategy", "adaptive");
    AddIpoptNumOption(nlp, "tol", 3.82e-6);
    AddIpoptStrOption(nlp, "mu_strategy", "adaptive");
    AddIpoptStrOption(nlp, "output_file", "ipopt.out");
    AddIpoptStrOption(nlp, "jacobian_approximation", "exact");
    AddIpoptStrOption(nlp, "hessian_approximation", "limited-memory");
    user_data.nlp = nlp;
    //SetIntermediateCallback(nlp, intermediate_cb);
    
    status = IpoptSolve(nlp, x, NULL, &obj, mult_g, mult_x_L, mult_x_U, &user_data);

Below there are the objective function along with the constraints:

static
bool eval_f(
   ipindex     n,
   ipnumber*   x,
   bool        new_x,
   ipnumber*   obj_value,
   UserDataPtr user_data
)
{
    
    fprintf(stderr, "Memory allocation failed.\n");
    struct MyUserData* optimization_data = user_data;
    double sum = 0.0;
    double* f = (double*) malloc((n) * sizeof(double));
    
    for(int i=0;i<n;i++){
        f[i]=(x[i+1]-x[i])/(optimization_data->ds[i]+(d(optimization_data->a,i)*d(x,i))+d(optimization_data->b,i));
    }
    for(int i=0;i<n;i++){
        sum+=pow(f[i],2)*optimization_data->ds[i];
    }
    *obj_value = sum;
    
    return true;
}
static
bool eval_g(
   ipindex     n,
   ipnumber*   x,
   bool        new_x,
   ipindex     m,
   ipnumber*   g,
   UserDataPtr user_data
)
{
   
    //g=(ipnumber*) malloc((m) * sizeof(ipnumber));
   
    struct MyUserData* optimization_data = user_data;
    double *kmap=(double*) malloc((n) * sizeof(double));
    double *velocity=(double*) malloc((n*sizeof(double)));
    double *map_h=(double*) malloc((n) * sizeof(double));
    double sum=0.0;
    
    double* f = (double*) malloc((n) * sizeof(double));
    for(int i=0;i<n;i++){
        f[i]=(x[i+1]-x[i])/(optimization_data->ds[i]+(d(optimization_data->a,i)*d(x,i))+d(optimization_data->b,i));
    }
    for(int i=0;i<n;i++){
        kmap[i]=x[i]*optimization_data->Kref;
        velocity[i]=sqrt(2*kmap[i]/optimization_data->m_map);
        map_h[i]=(optimization_data->fmapcoeff[0][0]*(exp(-pow((velocity[i]-optimization_data->fmapcoeff[1][0])/optimization_data->fmapcoeff[2][0],2)))+
        optimization_data->fmapcoeff[0][1]*(exp(-pow((velocity[i]-optimization_data->fmapcoeff[1][1])/optimization_data->fmapcoeff[2][1],2)))+
        optimization_data->fmapcoeff[0][2]*(exp(-pow((velocity[i]-optimization_data->fmapcoeff[1][2])/optimization_data->fmapcoeff[2][2],2))));
    }
    for (int i = 0; i < n; i++) {
        g[i]=f[i]-map_h[i];
    }
    for (int i = 0; i < n; i++) {
        sum += sqrt(1 / (2 * d(x, i))) * optimization_data->ds[i];
    }
    for(int i=0;i<n+1;i++){
        g[i]=sum;  
    }
    
    
    

    return true;
}

In addition there is the eval_grad_f function but i want to let ipopt to compute manually both the jacobians and hessians if possible.

static
bool eval_grad_f(
   ipindex     n,
   ipnumber*   x,
   bool        new_x,
   ipnumber*   grad_f,
   UserDataPtr user_data
)
{   
    double* aux1 = (double*)malloc(n * sizeof(double));
    double* aux2 = (double*)malloc(n * sizeof(double));
    struct MyUserData* optimization_data = user_data;
    double* f = (double*) malloc((n) * sizeof(double));
    
    for(int i=0;i<n;i++){
        f[i]=(x[i+1]-x[i])/(optimization_data->ds[i]+(d(optimization_data->a,i)*d(x,i))+d(optimization_data->b,i));
    }
    if (aux1 == NULL || aux2 == NULL) {
        fprintf(stderr, "Memory allocation failed.\n");
        exit(EXIT_FAILURE);
    }

    for (int i = 0; i < n; i++) {
        aux1[i] = 2 * optimization_data->ds[i] * f[i] * (0.5 * optimization_data->a[i] + 1.0 / optimization_data->ds[i]);
        aux2[i] = 2 * optimization_data->ds[i] * f[i] * (0.5 * optimization_data->a[i] - 1.0 / optimization_data->ds[i]);
    }

    for (int i = 0; i < n - 1; i++) {
        grad_f[i] = aux1[i] + aux2[i + 1];
    }
    grad_f[n - 1] = aux1[n - 1];

    // Free allocated memory
    free(aux1);
    free(aux2);
    return true;
}

Thank you in advance for your help and sorry for bothering you with this stuff :D, see ya! :)

[svigerske]

jacobian_approximation=exact means that the user has to provide the jacobian. Maybe there is an issue in your eval_jac_g.
The Jacobian and objective gradient approximations are not going to give decent performance, if they work at all.

In your eval_g, this doesn't look right:

    for(int i=0;i<n+1;i++){
        g[i]=sum;  
    }

Maybe just g[n] = sum;.

[Alfespo97]

Thanks for your answer, however that was a "debug" statement that i forgot to delete, i changed that line of code in g[n] = sum; but i still got segmentation fault. I have got also another question regarding the definition of constraints: i've got N non-linear inequality constraints such as f[i]-map_h[i]<=0 (i=0...N) and an equality constraint such as 1 / (2 * d(x, i))) * ds[i]==tmax (i=0...N) where x is the optimization variable, is the way i declared those constraints the right way or i forgot something? Could I share the definition of my model, instead of the code, to better understand the problem?

[svigerske]

Sounds right to me. f[i]-map_h[i]<=0 doesn't look non-linear to me.

A segmentation fault is usually due to an implementation error, looking at a model definition wouldn't help on that. Try some tool like valgrind for debugging.

[Alfespo97]

Ok i will try with valgrind to detect the error, thank you for your help, really appreciate it!