Main Content

最良実行可能点の取得

この例では、fmincon で見つかる最良実行可能点の取得方法を示します。

補助関数 bigtoleft は、x(1) 座標が負になるほど急速に負になる 3 次元変数 x 内の 3 次多項式目的関数です。その勾配は 3 要素のベクトルです。補助関数 bigtoleft のコードは、この例の終わりに掲載しています。

この例の制約集合は 2 つの円錐の内部交差部分です。1 つの円錐は上向きで、もう 1 つの円錐は下向きです。制約関数は各円錐に対して 1 つの要素を含む 2 要素のベクトルです。この例は 3 次元であるため、制約の勾配は 3 行 2 列の行列です。補助関数 twocone のコードは、この例の終わりに掲載しています。

目的関数を使用して関数 plottwoconecons を呼び出すことによって色分けされた制約の Figure を作成します (この例の最後にコードが掲載されています)。

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  
Objective function returned Inf; trying a new point...
   12         300  -5.691634e+258     1.355e+86     1.766e-43    1.538e+141    1.967e+173  

Feasible point with lower objective function value found, but optimality criteria not satisfied. See output.bestfeasible..


Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+02.

解および解法プロセスの検証

解、目的関数値、終了フラグおよび関数評価と反復の回数を調べます。

disp(x)
   1.0e+85 *

   -7.8891    7.7904   -2.4638
disp(fval)
 -5.6916e+258
disp(eflag)
     0
disp(output.constrviolation)
   1.3551e+86

目的関数値は、非常に大きな負の値である –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      -5.6916e+258         1.3551e+86     
    Best Feasible         -252.99                  0     

解の改善:範囲の設定

実行可能解を取得する方法は複数あります。1 つは変数で範囲を設定する方法です。

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, but optimality criteria not satisfied. See output.bestfeasible..


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 の最終値よりわずかに小さいです。最終点は、実行不可能性がより低く、1 次の最適性の尺度です。どちらの点がより適切でしょうか。これらはほぼ同じですが、それぞれに一方の点より優れた点があります。解の詳細を確認するために、表示形式を 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.000147102542200628 
    Best Feasible    -2894.12553454177    4.11274928779903e-07         0.21002299543909 

解の改善:別のアルゴリズムの使用

実行可能解を取得する別の方法は、境界のない '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, but optimality criteria not satisfied. See output.bestfeasible..


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 -3.500000000245461]
               fval: -2.894125047137579e+03
    constrviolation: 3.613823285064655e-08
      firstorderopt: 0.004069724066085

ここでも 2 つの解はほぼ同じですが、bestfeasible 解は終了前の反復で出力されます。最終解 xip のほうが 1 次の最適性の尺度と実行可能性が優れていますが、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

参考

関連するトピック