获取最佳可行点
此示例显示如何获取 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
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.
检查解和求解过程
检查解、目标函数值、退出标志以及函数计算和迭代的次数。
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
目标函数值小于-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 -5.6916e+258 1.3551e+86
Best Feasible -252.99 0
改进解：设置边界
有多种方法可以获得可行解。一种方法是设置变量的边界。
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.
迭代显示表明，fmincon 从可行点（可行性 0）开始，经过几次迭代不可行，再次达到 0 可行性，然后在剩余的迭代中具有较小但非零的不可行性。求解器再次报告，它在终点 xb 以外的点发现了较低的可行值。查看最终点和目标函数值，并报告目标函数值较低的可行点。
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
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.12500606454 4.61884486213648e-09 0.000147102542200628
Best Feasible -2894.12553454177 4.11274928779903e-07 0.21002299543909
改进解：使用另一种算法
获得可行解的另一种方法是使用 '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.
迭代显示再次表明 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('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
另请参阅
相关主题
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)