how to make fmincon honor constraints at each step?

Question

I am trying to estimate the parameters of a nonlinear filter using fmincon in the MATLAB Optimization Toolboox, keeping the eigenvalues of the state equation matrix less than one. The documentation for fmincon using the 'sqp' algorithm says that it respects the constraints at all iterations, and yet at the very first iteration of my code it tries a point implying unstable eigenvalues. I can't run the filter at these parameter values, and so my code crashes.

Does the documentation really mean that 'AlwaysHonorConstraints' only applies to the inequality constraints passed in as vectors (lb and ub)? Is there any way to ensure that the non-linear constraints are satisfied at each step?

Here is the relevant piece of my code:

c = @(x) max(abs(eig(reshape(x(8:16), 3, 3)))) - 1;
ceq = @(x) [];                 
nonlcon = @(x) deal(c(x), ceq(x));

obj = @(x) -unscented_kalman1(x, fulldiv, crsp_allEverything);

fminconOptions = optimset('Display', 'iter-detailed', 'Algorithm', 'sqp', ...
 'TolX', errTol, 'UseParallel', 'always', 'AlwaysHonorConstraints', 'bounds');

maxX = fmincon(obj, x0, [], [], [], [], [], [], nonlcon, fminconOptions);

Thanks!

Answer 1

The documentation indeed states that the 'AlwaysHonorConstraints' option applies to the bound constraints only, so lb and ub. In other words, fmincon considers the entire domain [lb ub] as the feasible region and ignores any and all linear and non-linear constraints.

Quite some time ago I wrote optimize, which is based on fminsearch and is available on the file exchange. Now, it's a long time ago, so I might misremember, but I remember implementing an option ('superstrict') that avoids any evaluation of the objective function outside the feasible region, including linear and non-linear constraints. You might want to give that a try.

But, sticking to fmincon, a possible workaround is to modify unscented_kalman1 such as to include penalties:

function unscented_kalman1(x, ...)

    %// evaluate constraint (yes, also here, unfortunately)
    c = @(x) max(abs(eig(reshape(x,2,2)))) - 1;
    penalties = c > 0;

    %// modify x by thresholding 
    if penalties          
        x = ... %// [MODIFY x SUCH THAT IT LIES INSIDE THE FEASIBLE REGION]
    end

    %// ... remainder of the function here

    %// then finally: 
    y = ... %// [THE OUTCOME OF THE FUNCTION];

    if penalties  
         y = y + c.^3;  %// ...or something similar
    end

end

This makes it possible to run your filter, but returning a much higher value for the objective function when the constraints are broken. This "forces" fmincon into the feasible region.

Take care though that the penalty function you use does not introduce discontinuities. Something like

 y = y + 1e8*c;

would not be suited, because if y is just inside the feasible region, the returned value is just y, but if it then just outside the region at the next iteration, it will suddenly jump up by millions. This is problematic, because SQP in fmincon uses first and second derivatives of y, which also turn to millions in the neighborhood of a discontinuity like that, screwing up the next iteration.