EXIT: Converged to a point of local infeasibility. Problem may be infeasible.

[wladimirpetrov]

Hello. Trying to solve the simple trajectory optimization task:

Find a set of N points on the 2D plane such that:
• the first and last point are positioned at (0, 0) and (10, 10), respectively;
• the distances between two consecutive points are equal;
• there is no point lying on a circle or radius 1 centered at (4, 4);
• the sum of the distances between consecutive points is minimized.

Here is the code:

#include "Ipopt/src/Interfaces/IpIpoptApplication.hpp"
#include "Ipopt/src/Interfaces/IpTNLP.hpp"
#include <cmath>
#include <vector>
#include <fstream>

using namespace Ipopt;

// Define a struct for 2D points
struct Point {
    double x;
    double y;
    Point(double x, double y) : x(x), y(y) {}
};
// predefined distance for the 1st iteration until the initial value isn't calculated
double distance = 1.0;

class PointSetProblem : public TNLP {
public:
    // Constructor
    PointSetProblem(int N, const Point& obstacle, const Point& pinit, const Point& pfin)
        : N_(N), obstacle_(obstacle), pinit_(pinit), pfin_(pfin) {}

    int getN() const {
        return N_;
    }

    // Define problem dimensions
    virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, Index& nnz_h_lag, IndexStyleEnum& index_style) {
        // Number of variables (2D coordinates for N points)
        n = 2 * N_;

        // Number of equality constraints (N-1 distances between points)
        m = 2 * N_ - 1;

        // Number of non-zero elements in the Jacobian of the constraints
        nnz_jac_g = 0;//2 * (N_ - 1);

        // Number of non-zero elements in the Hessian of the Lagrangian (0 if not available)
        nnz_h_lag = 0;

        // Use C-style indexing
        index_style = TNLP::C_STYLE;
        return true;
    }

    // Define bounds for variables and constraints
    virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u, Index m, Number* g_l, Number* g_u) {
    // Variable bounds (all points are in [0, 10] x [0, 10])
        for (int i = 0; i < N_; i++) {
            x_l[2 * i] = -1000;//-std::numeric_limits<double>::infinity(); // x coordinate (lower bound)
            x_u[2 * i] = 1000;//std::numeric_limits<double>::infinity(); // x coordinate (upper bound)
            x_l[2 * i + 1] = -1000;//-std::numeric_limits<double>::infinity(); // y coordinate (lower bound)
            x_u[2 * i + 1] = 1000;//std::numeric_limits<double>::infinity(); // y coordinate (upper bound)
            std::cout << "x_l[" << 2 * i << "] = " << x_l[2 * i] << std::endl;
            std::cout << "x_u[" << 2 * i << "] = " << x_u[2 * i] << std::endl;
            std::cout << "x_l[" << 2 * i + 1 << "] = " << x_l[2 * i + 1] << std::endl;
            std::cout << "x_u[" << 2 * i + 1 << "] = " << x_u[2 * i + 1] << std::endl;
        }

        // Additional constraints to enforce the initial and final points
        // Located exactly at (0, 0) and (10, 10)
        x_l[0] = x_u[0] = 0.0; // x coordinate of the initial point
        x_l[1] = x_u[1] = 0.0; // y coordinate of the initial point
        x_l[2 * (N_ - 1)] = x_u[2 * (N_ - 1)] = 10.0; // x coordinate of the final point
        x_l[2 * (N_ - 1) + 1] = x_u[2 * (N_ - 1) + 1] = 10.0; // y coordinate of the final point
        std::cout << "x_l[" << 0 << "] = " << x_l[0] << std::endl;
        std::cout << "x_u[" << 0 << "] = " << x_u[0] << std::endl;
        std::cout << "x_l[" << 1 << "] = " << x_l[1] << std::endl;
        std::cout << "x_u[" << 1 << "] = " << x_u[1] << std::endl;

        std::cout << "x_l[" << 2 * (N_ - 1) << "] = " << x_l[2 * (N_ - 1)] << std::endl;
        std::cout << "x_u[" << 2 * (N_ - 1) << "] = " << x_u[2 * (N_ - 1)] << std::endl;
        std::cout << "x_l[" << 2 * (N_ - 1) + 1 << "] = " << x_l[2 * (N_ - 1) + 1] << std::endl;
        std::cout << "x_u[" << 2 * (N_ - 1) + 1 << "] = " << x_u[2 * (N_ - 1) + 1] << std::endl;
        double radius = 1;

        // Constraint bounds (distances between points)
        for (int i = 0; i < 2*N_- 1; i++) {
            if (i == 0) {
                // For the first element of the equal distance i = 1
                g_l[i] = 0; // Lower bound for distance constraints
                g_u[i] = 1000;//std::numeric_limits<double>::infinity(); // Upper bound for distance constraints
                std::cout << "g_l[" << i << "] = " << g_l[i] << std::endl;
                std::cout << "g_u[" << i << "] = " << g_u[i] << std::endl;
            } else if (i >= 1 && i < N_ - 1) {
                // For the rest of the points for equal distance
                g_l[i] = 0; // Lower bound for distance constraints
                g_u[i] = 0; // Upper bound for distance constraints
                std::cout << "g_l[" << i << "] = " << g_l[i] << std::endl;
                std::cout << "g_u[" << i << "] = " << g_u[i] << std::endl;
            } else {
                // For Circle avoidence 
                g_l[i] = radius; // Lower bound for distance constraints
                g_u[i] = 1000;//std::numeric_limits<double>::infinity();
                //std::cout << "g_l = " << g_l << std::endl << std::flush;
                std::cout << "g_l[" << i << "] = " << g_l[i] << std::endl;
                std::cout << "g_u[" << i << "] = " << g_u[i] << std::endl; 

                
            }  
        
        }
        return true;
    }

    // Corrected initialization of the starting point
    virtual bool get_starting_point(Index n, bool init_x, Number* x, bool init_z, Number* z_L, Number* z_U, Index m, bool init_lambda, Number* lambda) {
        // Initialize points with equally spaced distances between pinit and pfin
        for (int i = 0; i < 2*N_; i++) {
            //x[i] = 0;
            x[0] = 0;
            x[1] = 0;
            x[2] = 2;
            x[3] = 2;
            x[4] = 4;
            x[5] = 4;
            x[6] = 6;
            x[7] = 6;
            x[8] = 8;
            x[9] = 8;
            x[10] = 10;
            x[11] = 10;
            
        }
        return true;
    }

    // Define the objective function
    virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value) {
        // Minimize the sum of distances between consecutive points
        obj_value = 0.0;
        for (int i = 0; i < N_ - 1; i++) {
            // Calculate the change in x and y between consecutive points
            double dx = x[2 * (i + 1)] - x[2 * i];
            double dy = x[2 * (i + 1) + 1] - x[2 * i + 1];
            // Calculate the Euclidean distance between consecutive points
            obj_value += std::sqrt(dx * dx + dy * dy);
        }
        return true;
    }

    // Define the constraints
    virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g) {
        // Calculate distances between consecutive points and use them as constraints
        for (int i = 0; i < N_- 1; i++) {
            if (i == 0) {
                // For i = 1
    
                double dx = x[2] - x[0];
                double dy = x[3] - x[1];
                g[i] = std::sqrt(dx * dx + dy * dy);
                std::cout << "g[" << i << "] = " << g[i] << std::endl;
            } else if (i >= 1 && i < N_ - 1) {
                std::cout << "x[" << i << "] = " << x[i] << std::endl;
                double dx_1 = x[2 * (i + 1)] - x[2 * i];
                std::cout << "dx_1[" << i << "] = " << dx_1 << std::endl;
                double dy_1 = x[2 * (i + 1)+1] - x[2 * i + 1];
                std::cout << "dy_1[" << i << "] = " << dy_1 << std::endl;
                double dx_0 = x[2] - x[0];
                std::cout << "dx_0[" << i << "] = " << dx_0 << std::endl;
                double dy_0 = x[3] - x[1];
                std::cout << "dy_0[" << i << "] = " << dy_0 << std::endl;
                g[i] = std::sqrt(dx_1 * dx_1 + dy_1 * dy_1)-std::sqrt(dx_0 * dx_0 + dy_0 * dy_0);
                std::cout << "g[" << i << "] = " << g[i] << std::endl;
                //std::cout << "x [" << i << "] = " << x[i] << std::endl;
            } 
        }
        
        for (int i = 0; i < N_; i++) {
            std::cout << "x[" << N_-1+i  << "] = " << x[N_-1+i ] << std::endl;
            std::cout << "obstacle_.x = " << obstacle_.x << std::endl;
            std::cout << "obstacle_.y = " << obstacle_.y << std::endl;
            double dx = x[2 * i] - obstacle_.x;
            std::cout << "dx[" << N_-1+i << "] = " << dx << std::endl;
            double dy = x[2 * i + 1] - obstacle_.y;
            std::cout << "dy[" << N_-1+i << "] = " << dy << std::endl;
            g[N_-1+i] = std::sqrt(dx * dx + dy * dy);
            std::cout << "g[" << N_-1+i << "] = " << g[N_-1+i] << std::endl;
            
        }
        return true;
    }

    // Define the Jacobian of the constraints
    virtual bool eval_jac_g(Index n, const Number* x, bool new_x, Index m, Index nele_jac, Index* iRow, Index* jCol, Number* values) {
        return true;
    }

    // Define the gradient of the objective function
    virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f) {
        return true;
    }

    // Method to finalize the solution
    virtual void finalize_solution(
        Ipopt::SolverReturn status, 
        Ipopt::Index n,
        const Ipopt::Number* x,
        const Ipopt::Number* z_L,
        const Ipopt::Number* z_U,
        Ipopt::Index m,
        const Ipopt::Number* g,
        const Ipopt::Number* lambda,
        Ipopt::Number obj_value,
        const Ipopt::IpoptData* ip_data,
        Ipopt::IpoptCalculatedQuantities* ip_cq

    ) {
        // Save the optimal points to a CSV file regardless of the optimization status
        std::ofstream outFile("optimal_points.csv");
        for (int i = 0; i < N_; i++) {
            outFile << x[2 * i] << "," << x[2 * i + 1] << std::endl;
        }
        outFile.close();
        std::cout << "Optimal points saved to 'optimal_points.csv'." << std::endl;

        if (status == Ipopt::SUCCESS) {
            std::cout << "Optimization succeeded." << std::endl;
        } else {
            std::cout << "IPOPT optimization status: " << status << std::endl;
        }
    }

private:
    int N_; // Number of points
    Point obstacle_; // Obstacle position
    Point pinit_; // Initial point
    Point pfin_; // Final point
    double distance;
};

int main() {
    int N = 6; // Number of points
    Point obstacle(4.0, 4.0);
    Point pinit(0.0, 0.0);
    Point pfin(10.0, 10.0);

    SmartPtr<TNLP> pointSetProblem = new PointSetProblem(N, obstacle, pinit, pfin);
    SmartPtr<IpoptApplication> app = IpoptApplicationFactory();

    // A smaller number pivots for sparsity, a larger number pivots for stability
    app->Options()->SetNumericValue("ma27_pivtol", 1e-8);//0.99 // 1e-8 // between 0 and 1
    // Ipopt may increase pivtol as high as ma27_pivtolmax to get a more accurate solution to the linear system
    app->Options()->SetNumericValue("ma27_pivtolmax", 0.0001);//0.99 // 0.0001 // between 0 and 1
    // The initial integer workspace memory = liw_init_factor * memory required by unfactored system. 
    // Ipopt will increase the workspace size by ma27_meminc_factor if required.
    app->Options()->SetNumericValue("ma27_liw_init_factor", 5.0); // 5.0 has to be 
    // The initial real workspace memory = la_init_factor * memory required by unfactored system. 
    // Ipopt will increase the workspace size by ma27_meminc_factor if required
    app->Options()->SetNumericValue("ma27_la_init_factor", 5.0); // 5.0
    // If the integer or real workspace is not large enough, Ipopt will increase its size by this factor.
    app->Options()->SetNumericValue("ma27_meminc_factor", 5.0); // 5.0

    app->Options()->SetStringValue("mu_strategy", "adaptive");
    app->Options()->SetStringValue("gradient_approximation", "finite-difference-values");

    // Set the Hessian approximation method to limited-memory
    app->Options()->SetStringValue("hessian_approximation", "limited-memory");

    // Adjust the maximum number of iterations
    app->Options()->SetIntegerValue("max_iter", 5000); // Change to your desired maximum iterations

    // Adjust the convergence tolerance
    app->Options()->SetNumericValue("tol", 1e-6); // Change to your desired tolerance

    app->RethrowNonIpoptException(true);
    ApplicationReturnStatus status = app->Initialize();
    if (status != Solve_Succeeded) {
        std::cout << "IPOPT initialization failed!" << std::endl;
        return -1;
    }

    status = app->OptimizeTNLP(pointSetProblem);

    if (status == Solve_Succeeded || status == Solved_To_Acceptable_Level) {
        Number obj_value;
        Number x[2 * N];
        pointSetProblem->eval_f(2 * N, x, true, obj_value);

        std::cout << "Optimal points:" << std::endl;
        for (int i = 0; i < N; i++) {
            std::cout << "Point " << i + 1 << ": (" << x[2 * i] << ", " << x[2 * i + 1] << ")" << std::endl;
        }
        std::cout << "Optimal Objective Value: " << obj_value << std::endl;
    } else {
        std::cout << "IPOPT optimization failed with status " << status << std::endl;
    }
    
    return 0;
}

The issue is that I am getting an error:

EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
IPOPT optimization status: 6
IPOPT optimization failed with status 2

The thing is that I'm pretty confident in the x and g bounds I assigned as well as I am confident in the x and g I provided. (this problem was initially solved by fmincon MATLAB function). As you can notice I print values and check the whether x and g satisfy the bounds. And as a matter of fact they do!

Maybe there are some issues with options of the IPOPT? May you help me in dealing with this problem?

[wladimirpetrov]

I tried to play with the

app->Options()->SetStringValue("mu_strategy", "monotone/adaptive"); (monotone mode led to more iterations 62 iterations in comparison to adaptive 15 iterations)
Tried to tune the options of the solver.
But nothing worked.

[svigerske]

Looking at your constraint, it doesn't seem correct that the Jacobian you provide is zero.

[svigerske]

Yes, you tell Ipopt that your Jacobian has 0 nonzero entries. So it assumes that all your constraint functions are constants. Therefore, Ipopt has no direction in which to reduce constraint violation, thus it declares local infeasibility.

If you want Ipopt to approximate the Jacobian, then you will have to enable the approximation of the Jacobian. In that case, you will still have to specify the nonzero pattern of the Jacobian, that is, specify a correct number for nnz_jac_g and implement eval_jac_g for the case it is called with values==NULL.

But, as said before, using approximated objective gradient and Jacobian is not recommended at all. If you cannot make use of some automatic differentiation code, then you may want to look at solvers for derivate-free optimization.

[wladimirpetrov]

Understood. However, Is there a way to optimize without specifying the Jacobian at all? Like I've done recently with much simpler task (with your help) link.
In that particular scenario the nnz_jac_g = 0 and I didn't provide Jacobian.

[svigerske]

#717 probably worked because your constraint was redundant.

If you cannot tell which variable appears in which constraint, then you can report a dense Jacobian to Ipopt (nnz_jac_g = n*m and fill out iCol and jRow accordingly), but you'll have to do this at least and enable the Jacobian approximation.

[wladimirpetrov]

Corrected the code by changing the line:
nnz_jac_g = n*m;
Added the following code to Jacobian part:

// Define the Jacobian of the constraints
        virtual bool eval_jac_g(Index n, const Number* x, bool new_x, Index m, Index nele_jac, Index* iRow, Index* jCol, Number* values) {
        if (values == NULL) {
            // Return the structure of the Jacobian by setting iRow and jCol
            for (Index i = 0; i < m; i++) {
                for (Index j = 0; j < n; j++) {
                    iRow[i * n + j] = i;
                    jCol[i * n + j] = j;
                }
            }
        }
        return true;
    }

Added the line:
app->Options()->SetStringValue("jacobian_approximation", "finite-differences");

EXIT: Maximum Number of Iterations Exceeded.

[wladimirpetrov]

may it be a problem that I am using original MA27 solver?

[svigerske]

Which problem? It doesn't say infeasible anymore. I wouldn't expect it to solve when using Jacobian or gradient approximations.