Получите лучшую допустимую точку

В этом примере показано, как получить лучшую допустимую точку, с которой сталкивается fmincon.

Функция помощника bigtoleft целевая функция кубического полинома в 3D переменной x это становится быстро отрицательным как x(1) координата становится отрицательной. Его градиент является трехэлементным вектором. Код для bigtoleft функция помощника появляется в конце этого примера.

Ограничительное множество для этого примера является пересечением внутренних частей двух конусов — один, обращенный вверх и один, обращенный вниз. Ограничительная функция является двухкомпонентным вектором, содержащим один компонент для каждого конуса. Поскольку этот пример 3D, градиент ограничения является 3-на-2 матрицей. Код для twocone функция помощника появляется в конце этого примера.

Создайте фигуру ограничений, окрашенных с помощью целевой функции путем вызова plottwoconecons функция, код которой появляется в конце этого примера.

figure1 = plottwoconecons;

Figure contains an axes object. The axes object contains 2 objects of type surface.

Создайте опции, чтобы использовать производные

Включить 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  -2.426511e+135     8.214e+44     1.000e+00     6.596e+44     1.167e+91  
   11          23  -7.428634e+134     4.785e+44     1.000e+00     2.543e+44     5.301e+90  
   12          25  -6.128293e+133     2.010e+44     1.000e+00     2.518e+44     1.003e+90  
   13          27  -5.589950e+131     9.249e+43     1.000e+00     1.509e+44     3.657e+88  
   14          34   2.149685e+130     1.006e+44     1.681e-01     5.456e+43     1.548e+88  
   15          36  -3.016774e+137     3.618e+45     1.000e+00     3.206e+45     2.907e+92  
   16          38  -3.628884e+137     3.577e+45     1.000e+00     4.039e+44     3.288e+92  
   17          87  -3.628884e+137     3.577e+45     2.569e-08     4.560e+37     3.288e+92  

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+45 *

   -3.3193   -0.0383    0.2577
disp(fval)
 -3.6289e+137
disp(eflag)
    -2
disp(output.constrviolation)
   3.5773e+45

Значение целевой функции меньше, чем –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      -3.6289e+137         3.5773e+45     
    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.


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.000147102542086941 
    Best Feasible    -2894.12553454177    4.11274928779903e-07        0.210022995438408 

Улучшите решение: используйте другой алгоритм

Другой способ получить возможное решение состоит в том, чтобы использовать '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: [-6.500000035892771 -7.634107877056255e-08 ... ]
               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

Смотрите также

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте