В этом примере показано, как получить лучшую допустимую точку, с которой сталкивается 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.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