Segmentation Fault of the code when choosing different value of the variables

[Lan-hhxx]

Hello everyone! I am trying to use the Idopt to solve the optimal question of the trajectory planning using C++ under the system Ubuntu20.04. And here is my code.
Code:

#include  <iostream>
#include <algorithm>
#include  <cassert>
#include  <cppad/ipopt/solve.hpp>
namespace {
	using CppAD::AD;
	class FG_eval {
	public:
		typedef CPPAD_TESTVECTOR( AD<double> ) ADvector;
		void operator()(ADvector& fg, const ADvector& x)
		{
			int N = 60;
			double T = 0.3;
			double g1 = 9.8;
			assert( fg.size() == 14 * (N - 1) + 1 );
			assert( x.size()  == 16* (N - 1) + 8 );
			fg[0] = 0;
			for (int i = 0; i < (N - 1); i++) {
				int m = 16 * i - 1;
                                int n = 16 * (i + 1) - 1;
				int j = 14 * i;
				fg[j + 1] = x[m + 9] - ((x[m + 1] + x[n + 1]) / 2 +  T * (x[m + 2] - x[n + 2]) / 8); 
         		fg[j + 2] = x[m + 10] - ((x[m + 2] + x[n + 2]) / 2 + T * (x[m + 7] - x[n + 7]) / 8); 
         		fg[j + 3] = x[m + 11] - ((x[m + 3] + x[n + 3]) / 2 + T * (x[m + 4] - x[n + 4]) / 8);
         		fg[j + 4] = x[m + 12] - ((x[m + 4] + x[n + 4]) / 2 + T * (x[m + 8] - x[n + 8]) / 8);
         		fg[j + 5] = x[m + 13] - ((x[m + 5] + x[n + 5]) / 2 + T * (x[m + 6] - x[n + 6]) / 8);
         		fg[j + 6] = x[m + 14] - ((x[m + 6] + x[n + 6]) / 2 + T * ((((-x[m+7]-2*x[m+4]*x[m+6]-g1*x[m+5])/x[m+3]) -(-x[n+7]-2*x[n+4]*x[n+6]-g1*x[n+5])/x[n+3])) / 8);
         		fg[j + 7] = x[m + 15] - ((x[m + 7] + x[n + 7]) / 2);
         		fg[j + 8] = x[m + 16] - ((x[m + 8] + x[n + 8]) / 2);
         		fg[j + 9] = x[m + 1] - x[n + 1] + T * (x[m + 2] + 4 * x[m + 10] + x[n + 2]) / 6;
         		fg[j + 10] = x[m + 2] - x[n + 2] + T * (x[m + 7] + 4 * x[m + 15] + x[n + 7]) / 6;
         		fg[j + 11] = x[m + 3] - x[n + 3] + T * (x[m + 4] + 4 * x[m + 12] + x[n + 4]) / 6;
         		fg[j + 12] = x[m + 4] - x[n + 4] + T * (x[m + 8] + 4 * x[m + 16] + x[n + 8]) / 6;
         		fg[j + 13] = x[m + 5] - x[n + 5] + T * (x[m + 6] + 4 * x[m + 14] + x[n + 6]) / 6;
				fg[j + 14] = x[m + 6] - x[n + 6] + T * ((-x[m+7]-2*x[m+4]*x[m+6]-g1*x[m+5])/x[m+3] + 4 * (-x[m+15]-2*x[m+12]*x[m+14]-g1*x[m+13])/x[m+11] + (-x[n+7]-2*x[n+4]*x[n+6]-g1*x[n+5])/x[n+3]) / 6;
			}
			return;
		}
	};
}

bool get_started(void)
{
	int N = 50;
	bool ok = true;
	size_t i;
	typedef CPPAD_TESTVECTOR( double ) Dvector;
	size_t nx = 16 * (N - 1) + 8;
	size_t ng = 14 * (N - 1);
	Dvector xi(nx);
	for(i = 0; i < 16 * (N - 1) + 8; i++)
	{
		xi[i] = 0.01;
	}
	double x0 = 0;double xf = 10;
	double l0 = 10;double lf = 0.5;
	double theta0 = 0;double thetaf = 0;
	double xv0 = 0;double xvf = 0;
	double lv0 = 0;double lvf = 0;
	double thetav0 = 0;double thetavf = 0;
	double xa0 = 0;double xaf = 0;
	double la0 = 0;double laf = 0;
	double x_min = std::min(x0, xf);double x_max = std::max(x0, xf);
	double l_min = std::min(l0, lf);double l_max = std::max(l0, lf);
	double theta_min = -0.5;double theta_max = 0.5;
	double xv_min = -1.1;double xv_max = 1.1;
	double lv_min = -2;double lv_max = 0;
	double thetav_min = -1;double thetav_max = 1;
	double xa_min = -0.5;double xa_max = 0.5;
	double la_min = -2;double la_max = 2;
	Dvector xl(nx), xu(nx);
	xl[0] = x0; xu[0] = x0;
	xl[1] = xv0; xu[1] = xv0;
	xl[2] = l0; xu[2] = l0;
	xl[3] = lv0; xu[3] = lv0;
	xl[4] = theta0; xu[4] = theta0;
	xl[5] = thetav0; xu[5] = thetav0;
	xl[6] = xa0; xu[6] = xa0;
	xl[7] = la0; xu[7] = la0;

	// nx = 16 * (N - `1`
	for(i = 0; i < (N - 1); i++)
	{
		int m = 16 * i - 1;
        int n = 16 * (i + 1) - 1;
		xl[m + 9] = x_min; xu[m + 9] = x_max;
		xl[m + 10] = xv_min; xu[m + 10] = xv_max;
		xl[m + 11] = l_min; xu[m + 11] = l_max;
		xl[m + 12] = lv_min; xu[m + 12] = lv_max;
		xl[m + 13] = theta_min; xu[m + 13] = theta_max;
		xl[m + 14] = thetav_min; xu[m + 14] = thetav_max;
		xl[m + 15] = xa_min; xu[m + 15] = xa_max;
		xl[m + 16] = la_min; xu[m + 16] = la_max;
		xl[n + 1] = x_min; xu[n + 1] = x_max;
		xl[n + 2] = xv_min; xu[n + 2] = xv_max;
		xl[n + 3] = l_min; xu[n + 3] = l_max;
		xl[n + 4] = lv_min; xu[n + 4] = lv_max;
		xl[n + 5] = theta_min; xu[n + 5] = theta_max;
		xl[n + 6] = thetav_min; xu[n + 6] = thetav_max;
		xl[n + 7] = xa_min; xu[n + 7] = xa_max;
		xl[n + 8] = la_min; xu[n + 8] = la_max;
	}
	
	int n = 16 * (N - 1);
	xl[n] = xf; xu[n] = xf;
	xl[n + 1] = xvf; xu[n + 1] = xvf;
	xl[n + 2] = lf; xu[n + 2] = lf;
	xl[n + 3] = lvf; xu[n + 3] = lvf;
	xl[n + 4] = thetaf; xu[n + 4] = thetaf;
	xl[n + 5] = thetavf; xu[n + 5] = thetavf;
	xl[n + 6] = xaf; xu[n + 6] = xaf;
	xl[n + 7] = laf; xu[n + 7] = laf;

	Dvector gl(ng), gu(ng);
	for(i = 0; i < 14 * (N - 1); i++)
	{
		gl[i] = 0; gu[i] = 0;
	}
	std::string options;
	// options += "Integer print_level  0\n";
	//  options += "String  sb           yes\n";
	// options += "Integer max_iter     5\n";
	// // approximate accuracy in first order necessary conditions;
	// // see Mathematical Programming, Volume 106, Number 1,
	// // Pages 25-57, Equation (6)
	//  options += "Numeric tol          1e-3\n";
	// options += "String  derivative_test            second-order\n";
	// maximum amount of random pertubation; e.g.,
	// when evaluation finite diff
	// options += "Numeric point_perturbation_radius  0.\n";
	CppAD::ipopt::solve_result<Dvector> solution;

	
	CppAD::ipopt::solve<Dvector, FG_eval>(
		options, xi, xl, xu, gl, gu, fg_eval, solution
	);
	ok &= solution.status == CppAD::ipopt::solve_result<Dvector>::success;
	for(i = 0; i < nx; i++)
	{	if (i % 16 == 1) //velocirty
        std::cout << "v[" << i << "] = " << solution.x[i] << std::endl;
	}
	return ok;
}

// main program that runs all the tests
int main(void)
{
    std::cout << "===== Ipopt with CppAD Testing =====" << std::endl;
    bool result = get_started();
    std::cout << "Final checking: " << result << std::endl;
}
// END C++

Result:
The solver will work and show the following information:

This is Ipopt version 3.14.14, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 532336
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 301476
Total number of variables............................: 776
variables with only lower bounds: 0
variables with lower and upper bounds: 776
variables with only upper bounds: 0
Total number of equality constraints.................: 686
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
.
.
.
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 0.0000000e+00 1.02e-02 3.68e+00 -1.7 1.86e+00 - 4.49e-01 2.02e-01f 1
61 0.0000000e+00 5.60e-03 2.27e+00 -1.7 1.30e+00 - 3.91e-01 5.19e-01f 1
62 0.0000000e+00 7.51e-03 1.38e-01 -1.7 8.30e-01 - 1.00e+00 1.00e+00f 1
63 0.0000000e+00 1.42e-05 1.53e-03 -2.5 3.78e-02 - 1.00e+00 1.00e+00h 1
64 0.0000000e+00 5.34e-07 1.92e-05 -3.8 8.91e-03 - 1.00e+00 1.00e+00h 1
65 0.0000000e+00 1.11e-07 1.25e-06 -5.7 1.60e-03 - 1.00e+00 1.00e+00h 1
66 0.0000000e+00 1.49e-09 1.87e-09 -8.6 1.84e-04 - 1.00e+00 1.00e+00h 1

Number of Iterations....: 66

                               (scaled)                 (unscaled)

Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 1.8730592004475365e-09 1.8730592004475365e-09
Constraint violation....: 6.8723325320688692e-10 1.4871727599397033e-09
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 3.4296130568355177e-09 3.4296130568355177e-09
Overall NLP error.......: 3.4296130568355177e-09 3.4296130568355177e-09

Number of objective function evaluations = 67
Number of objective gradient evaluations = 67
Number of equality constraint evaluations = 67
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 67
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 66
Total seconds in IPOPT = 136.355

EXIT: Optimal Solution Found.

And then I have encountered with the following question:
When I set the variable N under 64, like 60 or 61 or 62 , it will solve the problem as above and give the optimal result. However, once i set the N over 64, it will definitely show the "segmentation fault" without any other information of the fault and stop the program. I don't know the reason of it.Any ideas of how to get around this issue? Thanks.

[svigerske]

Since this is using CppAD's interface to Ipopt, you may want to try with the CppAD project.
Or just debug by yourself.

[Lan-hhxx]

ok thank you! I have built the ThirdParty-HSL and Ipopt during the installation. And I wonder where I should add the flag --enable-debug in my code to the call of configure. I am really not expert in debugging

[Lan-hhxx]

I have read the https://coin-or.github.io/Ipopt/OPTIONS.html but i still don't know how to change the option from ma27 to ma57 solver in my code. In the part std::string options; it seems to change the proceed parameter of the solution instead of the linear solver so i don't know where to switch the linear solver in my code.Could you please give some suggestion? And Merry Christmas!

[Lan-hhxx]

Besides, is there any limitation or option for the total number of variables of ipopt from the solver information cause i doubt once the total number is over 1000, it will lead to the segmentation fault.

[svigerske]

I have read the https://coin-or.github.io/Ipopt/OPTIONS.html but i still don't know how to change the option from ma27 to ma57 solver in my code. In the part std::string options; it seems to change the proceed parameter of the solution instead of the linear solver so i don't know where to switch the linear solver in my code.Could you please give some suggestion? And Merry Christmas!

I presume the way to do this with CppAD's Ipopt interface is

 options += "String  linear_solver ma57\n";

Besides, is there any limitation or option for the total number of variables of ipopt from the solver information cause i doubt once the total number is over 1000, it will lead to the segmentation fault.

By default, Ipopt's integer type, which is used to index variables, is a signed 32-bit integers. So there is some limit around 2^31, maybe less when slack variables are introduced.
Memory on the machine is usually a more limiting factor.

[Lan-hhxx]

Thank you! I have already solved this problem by adding the code options += "String linear_solver mumps\n" Thank you for your advice and instruction!