获取最佳可行点
此示例显示如何获取 fmincon 遇到的最佳可行点。
辅助函数 bigtoleft 是三维变量 x 中的三次多项式目标函数，随着 x(1) 坐标变为负值，该函数迅速变为负值。它的梯度是三元素向量。bigtoleft 辅助函数的代码出现在此示例末尾。
此示例的约束集是两个圆锥体（一个尖端向上，一个尖端向下）内部的交集。约束函数是一个双分量向量，每个锥包含一个分量。一个由于此示例是三维的，约束的梯度是一个 3×2 矩阵。twocone 辅助函数的代码出现在此示例末尾。
通过调用 plottwoconecons 函数创建使用目标函数着色的约束图形，该函数的代码出现在本示例的末尾。
figure1 = plottwoconecons;
创建使用导数的选项
为了使 fmincon 能够使用目标梯度和约束梯度，请设置适当的选项。选择 "sqp" 算法。
options = optimoptions("fmincon",Algorithm="sqp",... SpecifyConstraintGradient=true,SpecifyObjectiveGradient=true);
为了检查解过程，请设置选项以返回迭代显示。
options.Display = "iter";尽量减少使用衍生信息
设置初始点 x0 = [-1/20,-1/20,-1/20]。
x0 = -1/20*ones(1,3);
该问题没有线性约束或边界。将这些参量设置为 []。
A = []; b = []; Aeq = []; beq = []; lb = []; ub = [];
求解。
[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>
检查解和求解过程
检查解、目标函数值、退出标志以及函数计算和迭代的次数。
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
目标函数值小于-1e250，这是一个非常负的值。约束违反值大于 1e85，这是一个很大的数，但其量级仍然比目标函数值小得多。退出标志也表明返回的解是不可行的。
为了恢复 fmincon 遇到的最佳可行点及其目标函数值，显示 output.bestfeasible 结构体。
disp(output.bestfeasible)
x: [-2.9297 -0.1813 -0.1652]
fval: -252.9919
constrviolation: 0
firstorderopt: 258.5032
正如您在下一节中看到的，bestfeasible 点不是一个好的解。它只是 fmincon 在迭代中遇到的最佳可行点。在这种情况下，即使 bestfeasible 不是一个解，它也比返回的不可行解更好。
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
output.bestfeasible 结构包含 fmincon 在每次迭代结束时遇到的最佳可行点。在极少数情况下，fmincon 可能在某些中间计算中遇到一个可行点，其目标函数值甚至低于报告值。
改进解：设置界限
有多种方法可以获得可行解。一种方法是设置变量的边界。
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>
迭代显示表明，fmincon 从可行点（可行性 0）开始，经过几次迭代不可行，再次达到 0 可行性，然后在剩余的迭代中具有较小但非零的不可行性。求解器再次报告，它在终点 xb 以外的点发现了较低的可行值。查看最终点和目标函数值，并报告目标函数值较低的可行点。
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
bestfeasible 点处的约束违反值约为 4.113e-7。在迭代显示中，这种不可行性出现在迭代 11。该次迭代报告的目标函数值为 -2.894126e3，略小于最终值 -2.894125e3。最终点具有较低的不可行性和一阶最优性测度。哪个点比较好？它们几乎相同，但每一点都有其可取之处。要查看解的详细信息，请将显示格式设置为 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
改进解：使用另一种算法
获得可行解的另一种方法是使用 'interior-point' 算法，即使没有边界，
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>
迭代显示再次表明 fmincon 在寻找解时到达了不可行的点，并且 fmincon 再次发出一条消息，表示它遇到了具有较低目标函数值的可行点。
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
再次，这两个解几乎相同，并且 bestfeasible 解来自结束之前的迭代。最终解 xip 具有较好的一阶最优性测度和可行性，但 bestfeasible 解的目标函数值略低，且不不可行性不算太大。
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
最后，将格式重置为默认的 short。
format short辅助函数
以下代码创建 bigtoleft 辅助函数。
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
以下代码创建 twocone 辅助函数。
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
此代码创建绘制非线性约束的函数。
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
