静的解析のための for ループの作成

20 個の電子を導電体に置く問題を考えてみましょう (最適化変数を使用した静電学における制約付き非線形最適化を参照してください)。電子は、導電体内に存在する制約に従い、その総位置エネルギーを最小化するように整列します。この例は、本来入れ子形式の for ループとして表される位置エネルギーと、for ループを変換して静的解析を実行する方法に焦点を当てています。

静的解析のための修正なしの問題作成

N = 20 個の電子について、問題と最適化変数を作成します。

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;

導電体には以下の制約があります。

elecprob.Constraints.spherec = (x.^2 + y.^2 + (z+1).^2) <= 1;
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=\sum _{i<j}\frac{1}{\Vert electron(i)-electron(j)\Vert }.$$

目的関数を最適化式として定義します。

energy = optimexpr;
for ii = 1:(N-1)
    for jj = (ii+1):N
        energy = energy + ((x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2)^(-1/2);
    end
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 を呼び出して問題を解きます。求解時間を測定します。

tic
[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.

<stopping criteria details>

sol = struct with fields:
    x: [20×1 double]
    y: [20×1 double]
    z: [20×1 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.↵↵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.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 8.090884e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.'
            bestfeasible: [1×1 struct]
     objectivederivative: "reverse-AD"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

toc

Elapsed time is 14.483005 seconds.

求解時間は 10 秒以上です (結果は異なる場合があります)。

効率的な静的解析のための問題の修正

次のビデオで手順を示しています。

最適化式をより効率的に動作させるために、for ループを抽出してローカル関数に変換します。まず、energy = optimexpr; から最後の end までコードの行を選択します。次に、選択したコードの行を右クリックして [ローカル関数への変換] を選択します。

結果として得られる関数 (この例の終わりに掲載) を次のように編集します。

• 関数名を myenergy に変更します。

• コードの最初の行を次のように編集します。

energy = function myenergy(N, x, y, z);

• コードの 2 行目で、初期値を optimexpr から 0 に変更します。

静的解析を活用するには、fcn2optimexpr を呼び出すことで、エネルギーを目的関数として問題に配置します。

energy = fcn2optimexpr(@myenergy,N,x,y,z);
elecprob.Objective = energy;

solve を呼び出して問題を解きます。求解時間を測定します。

tic
[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.

<stopping criteria details>

sol = struct with fields:
    x: [20×1 double]
    y: [20×1 double]
    z: [20×1 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.↵↵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.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 8.090884e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.'
            bestfeasible: [1×1 struct]
     objectivederivative: "reverse-AD"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

toc

Elapsed time is 0.686725 seconds.

求解時間が大幅に短縮され、1 秒未満になります。

解析なしとの比較

目的関数をブラック ボックスとして扱うことで問題を解くこともできます。ブラック ボックスの目的関数は自動微分を活用できないので、ソルバーは有限差分ステップを使用して目的勾配を推定する必要があります。ただし、基となる関数の評価が高速な場合はブラック ボックスも高速に実行できます。fcn2optimexpr を設定してブラック ボックスの目的関数を作成した場合の求解時間を比較します。

energy = fcn2optimexpr(@myenergy,N,x,y,z,Analysis="off");
elecprob.Objective = energy;
tic
[sol,fval] = solve(elecprob,init)

Solving problem using fmincon.

Solver stopped prematurely.

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

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

fval = 163.3769

toc

Elapsed time is 0.392945 seconds.

関数の評価制限を超えるため、ソルバーが最適化を完了しません。また、結果の fval は前の 2 つの解より大きくなっています。関数の評価制限を増やし、再度実行します。

opts = optimoptions("fmincon",MaxFunctionEvaluations=1e4);
tic
[sol,fval,exitflag,output] = solve(elecprob,init,Options=opts)

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.

<stopping criteria details>

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

fval = 163.0099

exitflag = 
    OptimalSolution

output = struct with fields:
              iterations: 96
               funcCount: 5978
         constrviolation: 0
                stepsize: 5.5563e-06
               algorithm: 'interior-point'
           firstorderopt: 6.9362e-06
            cgiterations: 0
                 message: '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.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 6.902142e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.'
            bestfeasible: [1×1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

toc

Elapsed time is 0.603105 seconds.

求解時間は静的解析とほぼ同じで、目的関数の値も同様です。構造体 output は、ソルバーが有限差分ステップを使用して目的勾配を推定することを示しています。前のケースでは 200 回未満でしたが、ソルバーは何千回もの関数評価を行います。関数の評価に時間がかかる場合、このメソッドは非効率的である可能性があります。

結論

場合によっては、静的解析のための for ループを変換することで求解時間を節約できることがあります。for ループを関数に変換するのは非常に簡単です。ループを強調表示し、右クリックしてコンテキスト メニューを表示し、ローカル関数または関数ファイルを作成するオプションを選択し、新しい関数の細部を編集します。

どのような場合でも、fcn2optimexpr を使用して for ループの式を問題に配置します。これにより、静的解析の今後の改善を、コードを書き換えることなく自動的に取得できます。

補助関数

次のコードは補助関数 myenergy を作成します。

function energy = myenergy(N, x, y, z)
energy = 0;
for ii = 1:(N-1)
    for jj = (ii+1):N
        energy = energy + ((x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2)^(-1/2);
    end
end
end

参考

fcn2optimexpr

関連するトピック

• 静的解析のための for ループの制約の変換
• 問題ベースの最適化ワークフロー
• 方程式を解くための問題ベースのワークフロー