В этом примере показано, как получить лучшую допустимую точку, с которой сталкивается fmincon
.
Функция помощника bigtoleft
целевая функция кубического полинома в 3D переменной x
это становится быстро отрицательным как x(1)
координата становится отрицательной. Его градиент является трехэлементным вектором. Код для bigtoleft
функция помощника появляется в конце этого примера.
Ограничительное множество для этого примера является пересечением внутренних частей двух конусов — один, обращенный вверх и один, обращенный вниз. Ограничительная функция является двухкомпонентным вектором, содержащим один компонент для каждого конуса. Поскольку этот пример 3D, градиент ограничения является 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.89412
6e3
, который немного меньше окончательного значения -2.89412
5e3
. Конечная точка имеет более низкую недопустимость и меру оптимальности первого порядка. Какая точка лучше? Они - почти то же самое, но каждая точка имеет некоторое требование того, чтобы быть лучше. Чтобы видеть детали решения, установите формат отображения на 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