Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

静電学における制約付き非線形最適化、問題ベース

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

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

問題の幾何形状

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

z-|x|-|y|x2+y2+(z+1)21.

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

[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. The axes contains 2 objects of type surface.

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

問題の変数の定義

この問題では 20 個の電子を使用します。制約によって xy のそれぞれの値に –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;

上記の制約コマンドによって、10 個の制約のベクトルが作成されます。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 つの線形不等式として記述します。各制約コマンドは、20 個の制約のベクトルを返します。

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

目的関数の定義

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

energy=i<j1electron(i)-electron(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: '...'
            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. The axes 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.

関連するトピック