Obtain Best Feasible Point
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
12 25 -4.518887e+258 8.636e+85 1.000e+00 7.539e+85 1.766e+173
13 70 -4.518887e+258 8.636e+85 1.070e-07 3.756e+78 1.766e+173
Feasible point with lower objective function value found, but optimality criteria not satisfied. See output.bestfeasible..
Converged to an infeasible point.
fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.
<stopping criteria details>
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.6403 0.4014 -0.9857
disp(fval)
-4.5189e+258
disp(eflag)
-2
disp(output.constrviolation)
8.6365e+85
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 -4.5189e+258 8.6365e+85
Best Feasible -252.99 0
The output.bestfeasible structure contains the best feasible point that fmincon encountered at the end of each of its iterations. In rare cases, it is possible that fmincon encounters a feasible point in some intermediate calculation that has an even lower objective function value than it reports.
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.270030e+03 2.143e+00 7.000e-01 4.865e+00 2.804e+02
6 14 -2.538935e+03 2.855e-02 1.000e+00 2.229e+00 1.715e+03
7 16 -2.703318e+03 2.331e-02 1.000e+00 5.517e-01 2.521e+02
8 19 -2.845097e+03 0.000e+00 1.000e+00 1.752e+00 8.874e+01
9 21 -2.896933e+03 2.149e-03 1.000e+00 1.709e-01 1.608e+01
10 23 -2.894135e+03 7.953e-06 1.000e+00 1.039e-02 2.028e+00
11 25 -2.894126e+03 4.114e-07 1.000e+00 2.313e-03 2.100e-01
12 27 -2.894125e+03 4.618e-09 1.000e+00 2.450e-04 1.470e-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.
<stopping criteria details>
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.4498e-04 -3.5000]
fval: -2.8941e+03
constrviolation: 4.1139e-07
firstorderopt: 0.2100
The constraint violation at the bestfeasible point is about 4.113e-7. In the iterative display, this infeasibility 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.12500606352 4.61806859419767e-09 0.000146960594747725
Best Feasible -2894.12553469475 4.11390740140405e-07 0.21002648400258
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.
<stopping criteria details>
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 shortHelper 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("on"); % 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
