В этом примере показано, как получить лучшую допустимую точку, с которой сталкиваются fmincon.
Функция помощника bigtoleft является кубической полиномиальной целевой функцией в трехмерной переменной x который быстро растет отрицательным как x(1) координата становится отрицательной. Его градиент является трехэлементным вектором. Код для bigtoleft Функция helper появится в конце этого примера.
Ограничительное множество для этого примера является пересечение интерьеров двух конусов - одного , обращенного вверх, и одного , обращенного вниз. Функция ограничения является двухкомпонентным вектором, содержащим по одному компоненту для каждого конуса. Поскольку этот пример является трехмерным, градиент ограничения является матрицей 3 на 2. Код для twocone Функция helper появится в конце этого примера.
Создайте рисунок ограничений, окрашенных с помощью целевой функции путем вызова plottwoconecons function, чей код появляется в конце этого примера.
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
The 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