# Constrained function minimization with `fminsearch` and `fminlbfgs`

## Contents

- Constrained function minimization with `fminsearch` and `fminlbfgs`
- Unconstrained optimization
- Optimization with bound constraints
- Linear constraints
- Non-linear constraints
- Global optimization
- Different algorithm: `fminlbfgs`
- Supplying gradients to `fminlbfgs`
- See also

## Unconstrained optimization

first, define a test function:

        clc, rosen = @(x) ( (1-x(1))^2 + 105*(x(2)-x(1)^2)^2 ) /1e4;

This is the classical Rosenbrück function, which has a global minimum at . The function is relatively hard to minimize, because that minimum is located in a long narrow ``valley'':

        k = 0; range = -5:0.05:5;
        z = zeros(numel(range));
        for ii = range
            m = 0; k = k + 1;
            for jj = range
                m = m + 1;
                z(k, m) = rosen([ii, jj]);
            end
        end
        [y, x] = meshgrid(range, range);
        S = surf(x, y, z, 'linestyle', 'none'); view(-213, 38), axis tight

        shading interp, material metal, lighting gouraud, colormap('hot')
        light('style', 'local', 'position', [-3 0 5]);
        set(S, 'ambientstrength', 0.8)

Optimizing the fully unconstrained problem with `minimize` indeed finds the global minimum:

        solution = minimize(rosen, [3 3])

    solution =
        9.999779492786837e-001    9.999535409003786e-001

## Optimization with bound constraints

Imposing a lower bound on the variables gives

        [solution, fval] = minimize(rosen, [3 3], [],[], [],[], [2 2])

    solution =
        2.000000000528625e+000    3.999991976788512e+000
    fval =
        1.000000007819864e-004

in the figure, this looks like

        zz = z;   zz(x > 2 & y > 2) = inf;
        ZZ = z;   ZZ(x < 2 & y < 2) = inf;

        figure, hold on
        S(1) = surf(x, y, zz, 'linestyle', 'none', 'FaceAlpha', 0.2);
        S(2) = surf(x, y, ZZ, 'linestyle', 'none');
        plot3(solution(1), solution(2), fval+0.5, 'gx', ...
            'MarkerSize', 20,...
            'linewidth', 5)

        xlabel('X(1)'), ylabel('X(2)')
        view(-196, 38), grid on, axis tight

        shading interp, material metal, lighting gouraud, colormap('hot')
        light('style', 'local', 'position', [-3 0 5]);
        set(S, 'ambientstrength', 0.8);

Similarly, imposing an upper bound yields

        solution = minimize(rosen, [3 3], [],[], [],[], [],[0.5 0.5])

        zz = z;   zz(x < 0.5 & y < 0.5) = inf;
        ZZ = z;   ZZ(x > 0.5 & y > 0.5) = inf;

        figure, hold on
        S(1) = surf(x, y, zz, 'linestyle', 'none', 'FaceAlpha', 0.2);
        S(2) = surf(x, y, ZZ, 'linestyle', 'none');
        plot3(solution(1), solution(2), fval+0.5, 'gx', ...
            'MarkerSize', 20,...
            'LineWidth', 5);

        xlabel('X(1)'), ylabel('X(2)')
        view(201, 38), grid on, axis tight

        shading interp, material metal, lighting gouraud, colormap('hot')
        light('style', 'local', 'position', [-3 0 5]);
        set(S, 'ambientstrength', 0.8);

    solution =
        4.999999984348473e-001    2.499977662682127e-001

Minimize with fixed at 3. In this case, `minimize` simply removes the variable before `fminsearch` sees it, essentially reducing the dimensionality of the problem. This is particularly useful when the number of dimensions N becomes large.

        minimize(rosen, [3 3], [],[], [],[], [-inf 3], [inf 3])

    ans =
        1.731445312499997e+000    3.000000000000000e+000

## Linear constraints

You can use linear inequality or equality constraints. For example, with the constraints

,
,
with A = [+2 +1], b = -2 and Aeq = [+1 -1], beq = -2

`minimize()` finds the following result:

        [solution, fval] = minimize(rosen, [3;3], [2 1],-2, [1 -1],-2)

    solution =
       -1.333351464194111e+000
        6.666487473804754e-001
    fval =
        1.350896218141616e-002

These constraints look like the following:

        xinds = 2*x+y <= -2;
        zz = z;   zz( xinds ) = inf;
        ZZ = z;   ZZ(~xinds ) = inf;

        Ax = z;   Axinds = abs(2*x+y + 2) < 1e-3;
        [x1, sortinds] = sort(x(Axinds));
        Ax = Ax(Axinds); Ax = Ax(sortinds);
        y1 = y(Axinds); y1 = y1(sortinds);

        Aeq = z;   Aeqinds = abs(x-y + 2) < 1e-3;
        [x2, sortinds] = sort(x(Aeqinds));
        Aeq = Aeq(Aeqinds); Aeq = Aeq(sortinds);
        y2 = y(Aeqinds); y2 = y2(sortinds);

        figure, hold on
        l1 = line([x1(1:end-1)';x1(2:end)'],[y1(1:end-1)';y1(2:end)'],[Ax(1:end-1)';Ax(2:end)']);
        l2 = line([x2(1:end-1)';x2(2:end)'],[y2(1:end-1).';y2(2:end)'],[Aeq(1:end-1).';Aeq(2:end)']);
        S(1) = surf(x, y, zz, 'linestyle', 'none', 'FaceAlpha', 0.2);
        S(2) = surf(x, y, ZZ, 'linestyle', 'none');
        l3 = plot3(solution(1)+0.4, solution(2)+0.8, fval+0.5, 'gx',...
            'MarkerSize', 20,...
            'LineWidth', 5);

        set(l1, 'color', 'b', 'linewidth', 2);
        set(l2, 'color', 'k', 'linewidth', 2);

        view(150, 30), grid on, axis tight
        xlabel('X(1)', 'interpreter', 'LaTeX'); ylabel('X(2)', 'interpreter', 'LaTeX');
        k = legend([l1(1); l2(1); l3],'inequality $$A\mathbf{x} \leq -2$$', ...
            'equality $$A_{eq}\mathbf{x} = -2$$', 'Solution');
        set(k, 'interpreter', 'LaTeX', 'location', 'NorthWest');

        shading interp, material metal, lighting phong, colormap('autumn')
        light('style', 'local', 'position', [-3 0 5]);
        set(S, 'ambientstrength', 0.8);

## Non-linear constraints

Also general nonlinear constraints can be used. A simple example:

nonlinear inequality:

nonlinear equality :

        options = setoptimoptions(...
            'TolFun', 1e-6, ...
            'TolX'  , 1e-6, ...
            'MaxFunEvals', inf,...
            'MaxIter', 1e4);

        [sol, fval, exitflag, output] = minimize(rosen, [-3; 3], [],[], [],[],...
            [],[], @nonlcon, options);

Note that `nonlcon` is a subfunction, listed below.

These constraints look like the following:

        zz = z;   zz(sqrt(x.^2 + y.^2) <= 2)   = inf;
        ZZ = z;   ZZ(sqrt(x.^2 + y.^2) >= 2.2) = inf;
        zZ = z;   zZ(x.^2 + y.^3 >= 1.0 + 0.1) = inf;
                  zZ(x.^2 + y.^3 <= 1.0 - 0.1) = inf;

        xX = x(isfinite(zZ));  [xX, inds] = sort(xX);
        yY = y(isfinite(zZ));  yY = yY(inds);
        zZ = zZ(isfinite(zZ)); zZ = zZ(inds);

        figure, hold on
        S(1) = surf(x, y, zz, 'linestyle', 'none', 'FaceAlpha', 0.2);
        S(2) = surf(x, y, ZZ, 'linestyle', 'none');
        L = line([xX(1:end-1)';xX(2:end)'],[yY(1:end-1)';yY(2:end)'],[zZ(1:end-1)';zZ(2:end)']);
        l3 = plot3(sol(1)+0.4, sol(2)+0.5, fval+1, 'gx', 'MarkerSize', 20, 'linewidth', 5);

        set(L, 'linewidth', 2, 'color', 'b');
        view(150, 50), grid on, axis tight

        k = legend([S(2); L(1); l3],'non-linear inequality $$c(x) < 0$$', ...
            'non-linear equality $$c_{eq}(x) = 0$$', 'Solution');
        set(k, 'interpreter', 'LaTeX', 'location', 'NorthWest');

        shading interp, material metal, lighting phong, colormap('autumn')
        light('style', 'local', 'position', [-3 0 5]);
        set(S, 'ambientstrength', 0.8);

Note that the `output` structure contains a field `constrviolation`:

        output

    output =
             iterations: 2402
              algorithm: 'Nelder-Mead simplex direct search'
                message: [1x340 char]
           ObjfuncCount: 4425
        ConstrfuncCount: 4426
        constrviolation: [1x1 struct]

The contents of which shows that all constraints have been satisfied:

        output.constrviolation
        output.constrviolation.nonlin_eq{:}
        output.constrviolation.nonlin_ineq{:}

    ans =
          nonlin_eq: {2x1 cell}
        nonlin_ineq: {2x1 cell}
    ans =
         1
    ans =
        1.290913642648661e-008
    ans =
         0
    ans =
         0

## Global optimization

This is the 2D sine-envelope-sine function. It has a single global minimum at [0,0], where the function assumes a value of 0. As you can imagine, it is hard to find this minimum when the initial estimates is not very close to the minimum already:

        sinenvsin = @(x) 3*sum( (sin(sqrt(x(:).'*x(:))).^2 - 0.5)./(1 + 0.001*x(:).'*x(:)).^2 + 0.5, 1);

        figure, hold on
        k = 0; range = -10:0.1:10;
        z = zeros(numel(range));
        for ii = range
            m = 0; k = k + 1;
            for jj = range
                m = m + 1;
                z(k,m) = sinenvsin([ii jj]);
            end
        end
        [y, x] = meshgrid(range, range);
        S = surf(x, y, z, 'linestyle', 'none');

        axis equal, view(-148,24)

        shading interp, material shiny, lighting phong , colormap('autumn')
        light('style', 'local', 'position', [-3 0 5]);
        set(S, 'ambientstrength', 0.6);

`minimize()` provides rudimentary support for this type of problem. Omitting the initial value x0 will re-start `minimize()` several times at randomly chosen initial values in the interval [lb ub]:

        options = setoptimoptions(...
            'popsize', 1e2, 'maxfunevals', 1e4,  'maxiter', 1e2);

        [sol,fval] = minimize(sinenvsin, [], [],[], [],[], -[5 5], +[5 5], [],options)

    sol =
        8.819312347263519e-005    4.062807372839927e-005
    fval =
        2.831428808081071e-008

Naturally, these types of problems may also have constraints:

        [solution,fval] = minimize(sinenvsin, [], [2 1],-2, [1 -1],-2,...
            -[5;5], +[5;5], [], options);

        xinds = 2*x+y <= -2;
        zz = z;   zz( xinds ) = inf;
        ZZ = z;   ZZ(~xinds ) = inf;

        Ax = z;   Axinds = abs(2*x+y + 2) < 1e-3;
        [x1, sortinds] = sort(x(Axinds));
        Ax = Ax(Axinds); Ax = Ax(sortinds);
        y1 = y(Axinds); y1 = y1(sortinds);

        Aeq = z;   Aeqinds = abs(x-y + 2) < 1e-3;
        [x2, sortinds] = sort(x(Aeqinds));
        Aeq = Aeq(Aeqinds); Aeq = Aeq(sortinds);
        y2 = y(Aeqinds); y2 = y2(sortinds);

        figure, hold on
        S(1) = surf(x, y, zz, 'linestyle', 'none', 'FaceAlpha', 0.2);
        S(2) = surf(x, y, ZZ, 'linestyle', 'none');
        l1 = line([x1(1:end-1)';x1(2:end)'],[y1(1:end-1)';y1(2:end)'],[Ax(1:end-1)';Ax(2:end)']);
        l2 = line([x2(1:end-1)';x2(2:end)'],[y2(1:end-1).';y2(2:end)'],[Aeq(1:end-1).';Aeq(2:end)']);
        l3 = plot3(solution(1)+0.5, solution(2)+2.8, fval+2, 'gx', 'MarkerSize', 20, 'linewidth', 5);

        xlabel('X(1)', 'interpreter', 'LaTeX');
        ylabel('X(2)', 'interpreter', 'LaTeX');
        set(l1, 'color', 'r', 'linewidth', 2);
        set(l2, 'color', 'k', 'linewidth', 2);
        view(150, 30), grid on, axis tight
        k = legend([l1(1); l2(1); l3],'inequality $$A\mathbf{x} \leq -2$$', ...
            'equality $$A_{eq}\mathbf{x} = -2$$', 'Solution');
        set(k, 'interpreter', 'LaTeX', 'location', 'NorthWest');
        view(170,80);

        shading interp, material shiny, lighting phong, colormap('autumn')
        light('style', 'local', 'position', [-3 0 5]);
        set(S, 'ambientstrength', 0.6);

## Different algorithm: `fminlbfgs`

`minimize()` also supports `fminlbfgs`, a limited-memory, Broyden/Fletcher/Goldfarb/Shanno optimizer, implemented by Dirk-Jan Kroon. Unlike `fminsearch()`, this algorithm uses gradient information to improve the overall optimization speed:

        options = setoptimoptions(...
            'Algorithm', 'fminsearch');
        [solution, fval, exitflag, output] = minimize(...
            rosen, [-3 -18], [],[], [],[], [],[], [], options);

        solution, fval
        output.funcCount

        options = setoptimoptions(...
            'Algorithm', 'fminlbfgs');
        [solution, fval, exitflag, output] = minimize(...
            rosen, [-3 -18], [],[], [],[], [],[], [], options);

        solution, fval
        output.funcCount

    solution =
        1.000004532655050e+000    1.000006749626582e+000
    fval =
        5.836059191618274e-014
    ans =
       196
    solution =
        1.000012840436605e+000    1.000017577323999e+000
    fval =
        7.060245922779116e-013
    ans =
       162

As can be seen, `fminlbfgs()` indeed uses less funtion evaluations to come to a comparable answer.

The great advantage of this minimizer over `fminsearch()` is that `fminlbfgs()` can handle very large problems efficiently. For problems of higher dimensionality, `fminsearch()` has poor convergence properties compared to `fminlbfgs()`.

## Supplying gradients to `fminlbfgs`

`fminlbfgs()` needs gradient information of both the objective and non-linear constraint functions. `minimize()` estimate this information numerically via finite differences, but this can be costly for larger problems. Therefore, `minimize()` can also accept gradient information computed by the objective funcion:

        options = setoptimoptions(...
            'TolX', 1e-8,...
            'TolFun', 1e-8,...
            'FinDiffType', 'central',...
            'MaxFunEvals', inf,...
            'MaxIter', 1e4,...
            'GoalsExactAchieve', 0,...
            'Algorithm', 'fminlbfgs',...  % This option specifies that your
            'GradObj'  , 'on');           % objective function also provides
                                          % gradient information as its second
                                          % output argument

        [solution, fval, exitflag, output] = minimize(...
            @rosen_with_gradient, [-3 -3], [],[], [],[], [],[], @nonlcon, options);

        solution, fval
        output.ObjfuncCount
        output.ConstrfuncCount

In case of non-linearly constrained problems, also Jacobian information of the non-linear constraint function can be provided:

        options = setoptimoptions(...
            'TolX', 1e-10,...
            'TolFun', 1e-10,...
            'MaxFunEvals', inf,...
            'MaxIter', 1e4,...
            'GoalsExactAchieve', 0,...
            'Algorithm' , 'fminlbfgs',... % This option specifies that your
            'GradObj'   , 'on',...        % non-linear constraint function also
            'GradConstr', 'on');          % provides Jacobian information as its
                                          % third and fourth output arguments

        [solution, fval, exitflag, output] = minimize(...
            @rosen_with_gradient, [-3 -3], [],[], [],[], [],[], @nonlcon_with_Jacobian, options);

        solution, fval
        output.ObjfuncCount
        output.ConstrfuncCount

    solution =
        4.575590956156075e-001    1.338147855799337e+000
    fval =
        1.340811773195559e-002
    ans =
       508
    ans =
       509

## See also

setoptimoptions, fminlbfgs.

    end

    function [fVal, gVal] = rosen_with_gradient(x)
        fVal = ( (1-x(1))^2 + 105*(x(2)-x(1)^2)^2 ) /1e4;
        gVal = [ -2*(1-x(1)) - 4*105*x(1)*(x(2)-x(1)^2)
                 2*105*(x(2)-x(1)^2)
                 ]/1e4;
    end

    function [c, ceq] = nonlcon(x)
        c   = x(1)^2 + x(2)^2 - 2;
        ceq = 0.2*x(1)^2 + 0.4*x(2)^3 - 1;
    end

    function [c, ceq, c_Jac, ceq_Jac] = nonlcon_with_Jacobian(x)
        c   = x(1)^2 + x(2)^2 - 2;
        ceq = 0.2*x(1)^2 + 0.4*x(2)^3 - 1;

        c_Jac   = 2*x;
        ceq_Jac = [0.4*x(1); 1.2*x(2)^2];
    end

    solution =
        4.581537275933739e-001    1.337944382211403e+000
    fval =
        1.339032847543977e-002
    ans =
       420
    ans =
       421