最適化変数を使用した静電学における制約付き非線形最適化

ライブスクリプトを開く

20 個の電子を導電体に配置する静電学問題を考えてみましょう。電子は、導電体内に存在する制約に従い、その総位置エネルギーを最小化するように自身を配置します。すべての電子は最小値で導電体の境界上にあります。電子には区別がないため、この問題には固有の最小値がありません (1 つの解で電子を並べ替えれば、別の有効な解になります)。この例は、Dolan、Moré、Munson [1] の論文に基づきます。

最適化変数を使用する場合 (問題ベースのアプローチ)、この例の目的関数と非線形制約関数はすべて最適化変数および式でサポートされる演算です。そのため、solve は自動微分を使用して勾配を計算します。詳細は、Optimization Toolbox の自動微分を参照してください。自動微分を使用しない場合、この例はすぐに MaxFunctionEvaluations 許容誤差に到達して停止します。Symbolic Math Toolbox™ を使用した同等なソルバーベースの例については、Symbolic Math Toolbox を使用した勾配とヘッシアンの計算を参照してください。

問題の幾何形状

この例では、以下の不等式によって定義された導電体を使用します。座標が $(x, y, z)$ である各電子について、次のとおりとします。

$\begin{array}{l} z \leq - | x | - | y | \\ x^{2} + y^{2} + (z + 1)^{2} \leq 1 . \end{array}$

これらの制約は、球体の上の角錐があるような形状の導電体を形成します。導電体を表示するには、次のコードを入力します。

[X,Y] = meshgrid(-1:.01:1);
Z1 = -abs(X) - abs(Y);
Z2 = -1 - sqrt(1 - X.^2 - Y.^2);
Z2 = real(Z2);
W1 = Z1; W2 = Z2;
W1(Z1 < Z2) = nan; % Only plot points where Z1 > Z2
W2(Z1 < Z2) = nan; % Only plot points where Z1 > Z2
hand = figure; % Handle to the figure, for later use
set(gcf,'Color','w') % White background
surf(X,Y,W1,'LineStyle','none');
hold on
surf(X,Y,W2,'LineStyle','none');
view(-44,18)

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

形状の上と下の表面の間にわずかな隙間があります。この隙間は、形状を作成するために使用された一般的なプロットルーチンの人工物です。このルーチンは、他の面に触れる 1 面上の四角形の領域を消去します。

問題の変数の定義

この問題では 20 個の電子を使用します。制約によって $x$ と $y$ のそれぞれの値に –1 ～ 1 の範囲が指定され、 $z$ の値に –2 ～ 0 の範囲が指定されます。問題の変数を定義します。

N = 20;
x = optimvar('x',N,'LowerBound',-1,'UpperBound',1);
y = optimvar('y',N,'LowerBound',-1,'UpperBound',1);
z = optimvar('z',N,'LowerBound',-2,'UpperBound',0);
elecprob = optimproblem;

制約の定義

この問題には 2 種類の制約があります。1 つ目は、球面制約であり、各電子に対して別々に適用される単純な多項不等式です。この球面制約を定義します。

elecprob.Constraints.spherec = (x.^2 + y.^2 + (z+1).^2) <= 1;

上記の制約コマンドによって、N 個の制約のベクトルが作成されます。show を使用して制約ベクトルを表示します。

show(elecprob.Constraints.spherec)

  ((x.^2 + y.^2) + (z + 1).^2) <= arg_RHS

  where:

        arg2 = 1;
        arg1 = arg2(ones(1,20));
        arg_RHS = arg1(:);

問題の 2 つ目の制約のタイプは線形です。線形制約はさまざまな方法で表現できます。たとえば、関数 abs を使用して絶対値制約を表すことができます。この方法で制約を表現するには、MATLAB 関数を記述し、fcn2optimexpr を使用して式に変換します。詳細は、非線形関数から最適化式への変換を参照してください。微分可能関数のみを使用する望ましいアプローチでは、絶対値制約を 4 つの線形不等式として記述します。各制約コマンドは、N 個の制約のベクトルを返します。

elecprob.Constraints.plane1 = z <= -x-y;
elecprob.Constraints.plane2 = z <= -x+y;
elecprob.Constraints.plane3 = z <= x-y;
elecprob.Constraints.plane4 = z <= x+y;

目的関数の定義

目的関数はシステムの位置エネルギーです。これは、各電子ペアの距離の逆数の合計です。

$e n e r g y = \sum_{i < j} \frac{1}{‖ e l e c t r o n (i) - e l e c t r o n (j) ‖} .$

目的関数を最適化式として定義します。良好なパフォーマンスを得るために、目的関数をベクトル化された形式で記述します。詳細は、効率的な最適化問題の作成を参照してください。

energy = optimexpr(1);
for ii = 1:(N-1)
    jj = (ii+1):N;
    tempe = (x(ii) - x(jj)).^2 + (y(ii) - y(jj)).^2 + (z(ii) - z(jj)).^2;
    energy = energy + sum(tempe.^(-1/2));
end
elecprob.Objective = energy;

最適化の実行

[0,0,–1] でセンタリングされた半径 1/2 の球上で無作為に分布された電子で最適化を開始します。

rng default % For reproducibility
x0 = randn(N,3);
for ii=1:N
    x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;
    x0(ii,3) = x0(ii,3) - 1;
end
init.x = x0(:,1);
init.y = x0(:,2);
init.z = x0(:,3);

solve を呼び出して問題を解きます。

[sol,fval,exitflag,output] = solve(elecprob,init)

Solving problem using fmincon.

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.

sol = struct with fields:
    x: [20x1 double]
    y: [20x1 double]
    z: [20x1 double]

fval = 163.0099

exitflag = 
    OptimalSolution

output = struct with fields:
              iterations: 94
               funcCount: 150
         constrviolation: 0
                stepsize: 2.8395e-05
               algorithm: 'interior-point'
           firstorderopt: 8.1308e-06
            cgiterations: 0
                 message: 'Local minimum found that satisfies the constraints....'
            bestfeasible: [1x1 struct]
     objectivederivative: "reverse-AD"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

解の表示

解を導電体に点としてプロットします。

figure(hand)
plot3(sol.x,sol.y,sol.z,'r.','MarkerSize',25)
hold off

Figure contains an axes object. The axes object contains 3 objects of type surface, line.

電子は制約境界上に非常に均等に分布しています。多くの電子がエッジと角錐の頂点に位置します。

参考文献

[1] Dolan, Elizabeth D., Jorge J. Moré, and Todd S. Munson. "Benchmarking Optimization Software with COPS 3.0." Argonne National Laboratory Technical Report ANL/MCS-TM-273, February 2004.