This line below occurs six times in optim 1.6.2

nz = 200 * eps; # This is arbitrary. Accuracy of equality constraints.

I recommend that nz be controlled by the user through options.

To see the effect of nz, an example is provided.

fun=@(x) sum(x.^2); X0=[10 11]';
Aeq=[1 -1]; Beq=[0]; % X(1)-X(2)=0
Aieq=[]; Bieq=[];
lb=[1 1]'; ub=[12 12]';
[X, FVAL, CVG,  outp]= fmincon (fun, X0, Aieq, Bieq, Aeq, Beq, lb, ub, nonlcon=[]);

error: Initial parameters violate constraints.
error: called from
    _lm_feasible_ at line 93 column 5
    fmincon at line 456 column 24

% try to satisfy the constraint (it works this time but it is not always true)

fcn=@(x) Aeq*x-Beq;
[X0, FVAL, INFO, OUTPUT, FJAC] = fsolve(fcn, X0);
[X, FVAL, CVG,  outp]= fmincon (fun, X0, Aieq, Bieq, Aeq, Beq, lb, ub, nonlcon=[])

Depending on the complexity of fcn and options, fsolve may not result in a solution that will meet 200*eps (4e-16) in fmincon.  We can continue to reduce TolX & TolFun for fdolve.  But a possible approach seems to be to relax nz or allow users to adjust nz as they wish (like TolX & TolFun).  I recommend the latter.  Thanks.