В этом примере показано, как получить наиболее возможную точку, встречающуюся при 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.
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.
Проверьте решение, значение целевой функции, флаг выхода и количество оценок и итераций функции.
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
Существует несколько способов получения выполнимого решения. Один из способов - установить границы для переменных.
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.538947e+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.
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.12500606462 4.61890525826902e-09 0.000147112019817541
Best Feasible -2894.1255345325 4.11267924604886e-07 0.21002297099706
Другой способ получить выполнимое решение - использовать '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.
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: [1x3 double]
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