Main Content

Obtain Best Feasible Point

Open Live Script

This example shows how to obtain the best feasible point encountered by fmincon.

The helper function bigtoleft is a cubic polynomial objective function in a three-dimensional variable x that grows rapidly negative as the x(1) coordinate becomes negative. Its gradient is a three-element vector. The code for the bigtoleft helper function appears at the end of this example.

The constraint set for this example is the intersection of the interiors of two cones—one pointing up, and one pointing down. The constraint function is a two-component vector containing one component for each cone. Because this example is three-dimensional, the gradient of the constraint is a 3-by-2 matrix. The code for the twocone helper function appears at the end of this example.

Create a figure of the constraints colored using the objective function by calling the plottwoconecons function, whose code appears at the end of this example.

figure1 = plottwoconecons;

Create Options To Use Derivatives

To enable fmincon to use the objective gradient and constraint gradients, set appropriate options. Choose the 'sqp' algorithm.

options = optimoptions('fmincon','Algorithm','sqp',...
    "SpecifyConstraintGradient",true,"SpecifyObjectiveGradient",true);

To examine the solution process, set options to return iterative display.

options.Display = 'iter';

Minimize Using Derivative Information

Set the initial point x0 = [-1/20,-1/20,-1/20].

x0 = -1/20*ones(1,3);

The problem has no linear constraints or bounds. Set those arguments to [].

A = [];
b = [];
Aeq = [];
beq = [];
lb = [];
ub = [];

Solve the problem.

[x,fval,eflag,output] = fmincon(@bigtoleft,x0,...
           A,b,Aeq,beq,lb,ub,@twocone,options);
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           1   -1.625000e-03     0.000e+00     1.000e+00     0.000e+00     8.250e-02  
    1           3   -2.490263e-02     0.000e+00     1.000e+00     8.325e-02     5.449e-01  
    2           5   -2.529919e+02     0.000e+00     1.000e+00     2.802e+00     2.585e+02  
    3           7   -6.408576e+03     9.472e+00     1.000e+00     1.538e+01     1.771e+03  
    4           9   -1.743599e+06     5.301e+01     1.000e+00     5.991e+01     9.216e+04  
    5          11   -5.552305e+09     1.893e+03     1.000e+00     1.900e+03     1.761e+07  
    6          13   -1.462524e+15     5.632e+04     1.000e+00     5.636e+04     8.284e+10  
    7          15   -2.573346e+24     1.471e+08     1.000e+00     1.471e+08     1.058e+17  
    8          17   -1.467510e+41     2.617e+13     1.000e+00     2.617e+13     1.789e+28  
    9          19   -8.716877e+72     2.210e+24     1.000e+00     2.210e+24     2.387e+49  
   10          21  -5.598248e+134     4.090e+44     1.000e+00     4.090e+44     4.368e+90  
   11          23  -5.691634e+258     1.355e+86     1.000e+00     1.136e+86    1.967e+173  
Objective function returned Inf; trying a new point...
   12         300  -5.691634e+258     1.355e+86     1.766e-43    1.538e+141    1.967e+173  

Feasible point with lower objective function value found, but optimality criteria not satisfied. See output.bestfeasible..


Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+02.

Examine Solution and Solution Process

Examine the solution, objective function value, exit flag, and number of function evaluations and iterations.

disp(x)
   1.0e+85 *

   -7.8891    7.7904   -2.4638

disp(fval)
 -5.6916e+258

disp(eflag)
     0

disp(output.constrviolation)
   1.3551e+86

The objective function value is smaller than –1e250, a very negative value. The constraint violation is larger than 1e85, a large amount that is still much smaller in magnitude than the objective function value. The exit flag also shows that the returned solution is infeasible.

To recover the best feasible point that fmincon encounters, along with its objective function value, display the output.bestfeasible structure.

disp(output.bestfeasible)
                  x: [-2.9297 -0.1813 -0.1652]
               fval: -252.9919
    constrviolation: 0
      firstorderopt: 258.5032

The bestfeasible point is not a good solution, as you see in the next section. It is simply the best feasible point that fmincon encountered in its iterations. In this case, even though bestfeasible is not a solution, it is a better point than the returned infeasible solution.

table([fval;output.bestfeasible.fval],...
    [output.constrviolation;output.bestfeasible.constrviolation],...
    'VariableNames',["Fval" "Constraint Violation"],'RowNames',["Final Point" "Best Feasible"])
ans=2×2 table
                         Fval        Constraint Violation
                     ____________    ____________________

    Final Point      -5.6916e+258         1.3551e+86     
    Best Feasible         -252.99                  0

Improve Solution: Set Bounds

There are several ways to obtain a feasible solution. One way is to set bounds on the variables.

lb = -10*ones(3,1);
ub = -lb;
[xb,fvalb,eflagb,outputb] = fmincon(@bigtoleft,x0,...
           A,b,Aeq,beq,lb,ub,@twocone,options);
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           1   -1.625000e-03     0.000e+00     1.000e+00     0.000e+00     8.250e-02  
    1           3   -2.490263e-02     0.000e+00     1.000e+00     8.325e-02     5.449e-01  
    2           5   -2.529919e+02     0.000e+00     1.000e+00     2.802e+00     2.585e+02  
    3           7   -4.867942e+03     5.782e+00     1.000e+00     1.151e+01     1.525e+03  
    4           9   -1.035980e+04     3.536e+00     1.000e+00     9.587e+00     1.387e+03  
    5          12   -5.270039e+03     2.143e+00     7.000e-01     4.865e+00     2.804e+02  
    6          14   -2.538946e+03     2.855e-02     1.000e+00     2.229e+00     1.715e+03  
    7          16   -2.703320e+03     2.330e-02     1.000e+00     5.517e-01     2.521e+02  
    8          19   -2.845099e+03     0.000e+00     1.000e+00     1.752e+00     8.873e+01  
    9          21   -2.896934e+03     2.150e-03     1.000e+00     1.709e-01     1.608e+01  
   10          23   -2.894135e+03     7.954e-06     1.000e+00     1.039e-02     2.028e+00  
   11          25   -2.894126e+03     4.113e-07     1.000e+00     2.312e-03     2.100e-01  
   12          27   -2.894125e+03     4.619e-09     1.000e+00     2.450e-04     1.471e-04  

Feasible point with lower objective function value found, but optimality criteria not satisfied. See output.bestfeasible..


Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

The iterative display shows that fmincon starts at a feasible point (feasibility 0), spends a few iterations infeasible, again reaches 0 feasibility, then has small but nonzero infeasibility for the remaining iterations. The solver again reports that it found a lower feasible value at a point other than the final point xb. View the final point and objective function value, and the reported feasible point with lower objective function value.

disp(xb)
   -6.5000   -0.0000   -3.5000

disp(fvalb)
  -2.8941e+03

disp(outputb.bestfeasible)
                  x: [-6.5000 2.4500e-04 -3.5000]
               fval: -2.8941e+03
    constrviolation: 4.1127e-07
      firstorderopt: 0.2100

The constraint violation at the bestfeasible point is about 4.113e-7. In the iterative display, this infeasibiliity occurs at iteration 11. The reported objective function value at that iteration is -2.894126e3, which is slightly less than the final value of -2.894125e3. The final point has lower infeasibility and first-order optimality measure. Which point is better? They are nearly the same, but each point has some claim to being better. To see the solution details, set the display format to long.

format long
table([fvalb;outputb.bestfeasible.fval],...
    [outputb.constrviolation;outputb.bestfeasible.constrviolation],...
    [outputb.firstorderopt;outputb.bestfeasible.firstorderopt],...
    'VariableNames',["Function Value" "Constraint Violation" "First-Order Optimality"],...
    'RowNames',["Final Point" "Best Feasible"])
ans=2×3 table
                      Function Value      Constraint Violation    First-Order Optimality
                     _________________    ____________________    ______________________

    Final Point      -2894.12500606454    4.61884486213648e-09     0.000147102542200628 
    Best Feasible    -2894.12553454177    4.11274928779903e-07         0.21002299543909

Improve Solution: Use Another Algorithm

Another way to obtain a feasible solution is to use the 'interior-point' algorithm, even without bounds,

lb = [];
ub = [];
options.Algorithm = "interior-point";
[xip,fvalip,eflagip,outputip] = fmincon(@bigtoleft,x0,...
           A,b,Aeq,beq,lb,ub,@twocone,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1   -1.625000e-03    0.000e+00    7.807e-02
    1       2   -2.374253e-02    0.000e+00    5.222e-01    8.101e-02
    2       3   -2.232989e+02    0.000e+00    2.379e+02    2.684e+00
    3       4   -3.838433e+02    1.768e-01    3.198e+02    5.573e-01
    4       5   -3.115565e+03    1.810e-01    1.028e+03    4.660e+00
    5       6   -3.143463e+03    2.013e-01    8.965e+01    5.734e-01
    6       7   -2.917730e+03    1.795e-02    6.140e+01    5.231e-01
    7       8   -2.894095e+03    0.000e+00    9.206e+00    1.821e-02
    8       9   -2.894107e+03    0.000e+00    2.500e+00    3.899e-03
    9      10   -2.894142e+03    1.299e-05    2.136e-03    1.371e-02
   10      11   -2.894125e+03    3.614e-08    4.070e-03    1.739e-05
   11      12   -2.894125e+03    0.000e+00    5.994e-06    5.832e-08

Feasible point with lower objective function value found, but optimality criteria not satisfied. See output.bestfeasible..


Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

The iterative display again shows fmincon reaching points that are infeasible in its search for a solution, and fmincon again issues a message that it encountered a feasible point with lower objective function value.

       disp(xip)
  -6.499999996950366  -0.000000032933162  -3.500000000098132

disp(fvalip)
    -2.894124995999976e+03

disp(outputip.bestfeasible)
                  x: [-6.500000035892771 -7.634107877056255e-08 -3.500000000245461]
               fval: -2.894125047137579e+03
    constrviolation: 3.613823285064655e-08
      firstorderopt: 0.004069724066085

Again, the two solutions are nearly the same, and the bestfeasible solution comes from the iteration before the end. The final solution xip has better first-order optimality measure and feasibility, but the bestfeasible solution has slightly lower objective function value and an infeasibility that it not too large.

table([fvalip;outputip.bestfeasible.fval],...
    [outputip.constrviolation;outputip.bestfeasible.constrviolation],...
    [outputip.firstorderopt;outputip.bestfeasible.firstorderopt],...
    'VariableNames',["Function Value" "Constraint Violation" "First-Order Optimality"],...
    'RowNames',["Final Point" "Best Feasible"])
ans=2×3 table
                      Function Value      Constraint Violation    First-Order Optimality
                     _________________    ____________________    ______________________

    Final Point      -2894.12499599998                       0     5.99383553128062e-06 
    Best Feasible    -2894.12504713758    3.61382328506465e-08      0.00406972406608475

Finally, reset the format to the default short.

format short

Helper Functions

This code creates the bigtoleft helper function.

function [f gradf] = bigtoleft(x)
% This is a simple function that grows rapidly negative
% as x(1) becomes negative
%
f = 10*x(:,1).^3 + x(:,1).*x(:,2).^2 + x(:,3).*(x(:,1).^2 + x(:,2).^2);

if nargout > 1

   gradf=[30*x(1)^2+x(2)^2+2*x(3)*x(1);
       2*x(1)*x(2)+2*x(3)*x(2);
       (x(1)^2+x(2)^2)];

end
end

This code creates the twocone helper function.

function [c ceq gradc gradceq] = twocone(x)
% This constraint is two cones, z > -10 + r
% and z < 3 - r

ceq = [];
r = sqrt(x(1)^2 + x(2)^2);
c = [-10+r-x(3);
    x(3)-3+r];

if nargout > 2

    gradceq = [];
    gradc = [x(1)/r,x(1)/r;
       x(2)/r,x(2)/r;
       -1,1];

end
end

This code creates the function that plots the nonlinear constraints.

function figure1 = plottwoconecons
% Create figure
figure1 = figure;

% Create axes
axes1 = axes('Parent',figure1);
view([-63.5 18]);
grid('on');
hold('all');

% Set up polar coordinates and two cones
r = linspace(0,6.5,14);
th = 2*pi*linspace(0,1,40);
x = r'*cos(th);
y = r'*sin(th);
z = -10+sqrt(x.^2+y.^2);
zz = 3-sqrt(x.^2+y.^2);

% Evaluate objective function on cone surfaces
newxf = reshape(bigtoleft([x(:),y(:),z(:)]),14,40)/3000;
newxg = reshape(bigtoleft([x(:),y(:),z(:)]),14,40)/3000;

% Create lower surf with color set by objective
surf(x,y,z,newxf,'Parent',axes1,'EdgeAlpha',0.25);

% Create upper surf with color set by objective
surf(x,y,zz,newxg,'Parent',axes1,'EdgeAlpha',0.25);
axis equal
xlabel 'x(1)'
ylabel 'x(2)'
zlabel 'x(3)'
end

See Also

Related Topics