Roundoff error when adding an equality constraint to SLSQP

[JeremyBois]

Hi,

Description

I'm trying to use the SLSQP optimizer with equality and inequality constraints to port a code from python (using scipy) to C# (using NLoptNet nlopt wrapper).

When using SLSQP version of scipy I can add equality constraints without any issues. When using SLSQP version from nlopt I get a NLOPT_ROUNDOFF_LIMITED error code.

Additional information

I search some information (in nlopt documentation) about the SLSQP optimizer implementation and both scipy and nlopt version seem to be based on the same fortran code.

I'm using a C# wrapper in my code (see issue BrannonKing/NLoptNet#16) but the same issue is raised with C++ version (see below).

Maybe related with #350

Minimal working code

#include <stdio.h>
#include <math.h>
#include <vector>

#include "nlopt.hpp"


int count = 0;

double myfunc(unsigned n, const double* variables, double* gradient, void* my_func_data)
{
    const float h = 0.001;
    auto f0 = pow(variables[1] * variables[0], 2);
    auto fi0 = pow(variables[1] * (variables[0] + h), 2);
    auto fi1 = pow((variables[1] + h) * variables[0], 2);
    if (gradient)
    {
        gradient[0] = (fi0 - f0) / h;
        gradient[1] = (fi1 - f0) / h;

        printf("> Iteration #%d: x = [%f, %f], gradient = [%f, %f]\n", count, variables[0], variables[1], gradient[0], gradient[1]);
    }
    count++;

    return f0;
}

double EqConstraint(unsigned int n, const double* variables, double* gradient, void* data)
{
    return 0;
}
double IEqConstraint(unsigned int n, const double* variables, double* gradient, void* data)
{
    return 0;
}


int main()
{
    nlopt::opt opt(nlopt::LD_SLSQP, 2);
    opt.set_maxeval(20);
    opt.set_xtol_rel(1e-4);
    opt.set_min_objective(myfunc, NULL);

    std::vector<double> lower_bounds(2);
    lower_bounds[0] = 1.0;
    lower_bounds[1] = 1.0;
    opt.set_lower_bounds(lower_bounds);

    std::vector<double> upper_bounds(2);
    upper_bounds[0] = 10.0;
    upper_bounds[1] = 10.0;
    opt.set_upper_bounds(upper_bounds);

    // Constraints
    double zero = 0.0;

    // @BUG Adding equality constraint breaks the code 
    opt.add_equality_constraint(EqConstraint, &zero, 0.001);
    // @BUG Adding equality constraint breaks the code 

    // Adding inequality constraint works
    opt.add_inequality_constraint(IEqConstraint, &zero, 0.001);
    
    try
    {
        // Let's go !
        double minf;
        std::vector<double> initialValues(2);
        initialValues[0] = 1.234; initialValues[1] = 5.678;

        auto result = opt.optimize(initialValues, minf);
    }
    catch (const std::exception&)
    {

    }
    
    printf("result: %d\n", opt.last_optimize_result());
}