ODE を並列に最適化する

このページは機械翻訳を使用して翻訳されました。元の英語を参照するには、ここをクリックします。

この例では、ODE のパラメータを最適化する方法を示します。

また、ODE ソリューションが目的関数と非線形制約関数の両方を返す場合に、それらの関数を 2 回計算しないようにする方法も示します。この例では、ソルバーの実行時間とソリューションの品質の観点から、patternsearch と ga を比較します。

並列コンピューティングを使用するには、Parallel Computing Toolbox™ ライセンスが必要です。

手順 1: 問題を定義します。
手順 2: ODE モデルを定式化します。
手順 3: パターン検索を使用して解決します。
手順 4: コストのかかるサブルーチンを 2 回呼び出すことは避けてください。
手順 5: 並列で計算します。
手順 6.遺伝的アルゴリズムと比較してください。

手順 1: 問題を定義します。

問題は、大砲の位置と角度を変えて、できるだけ壁を越えて発射物を発射することです。大砲の砲口速度は300m/sです。壁の高さは20メートルです。大砲が壁に近すぎると、急な角度で発射しなければならなくなり、発射物が十分遠くまで飛びません。大砲が壁から離れすぎている場合も、発射物は十分な距離を飛ばないことになります。

空気抵抗により発射物の速度が遅くなります。抵抗力は速度の二乗に比例し、比例定数は 0.02 です。重力は発射体に作用し、一定の 9.81 m/s² で下向きに加速します。したがって、軌道x(t)の運動方程式は

$\frac{d^{2} x (t)}{d t^{2}} = - 0.02 ‖ v (t) ‖ v (t) - (0, 9.81),$

ここで、 $v (t) = d x (t) / d t .$ です。

初期位置 x0 と初期速度 xp0 は 2 次元ベクトルです。ただし、初期の高さ x0(2) は 0 なので、初期位置はスカラー x0(1) のみによって決まります。また、初期速度 xp0 は大きさが 300 (銃口速度) なので、初期角度 (スカラー) のみに依存します。初期角度 th の場合、 xp0 = 300*(cos(th),sin(th)) です。したがって、最適化問題は 2 つのスカラーのみに依存するため、2 次元の問題になります。水平距離と角度を決定変数として使用します。

手順 2: ODE モデルを定式化します。

ODE ソルバーでは、モデルを 1 次システムとして定式化する必要があります。軌道ベクトル (x₁(t),x₂(t)) にその時間微分 (x'₁(t),x'₂(t)) を加えて、4 次元軌道ベクトルを形成します。この拡大ベクトルに関して、微分方程式は次のようになる。

$\frac{d}{d t} x (t) = [\begin{matrix} x_{3} (t) \\ x_{4} (t) \\ - .02 ‖ (x_{3} (t), x_{4} (t)) ‖ x_{3} (t) \\ - .02 ‖ (x_{3} (t), x_{4} (t)) ‖ x_{4} (t) - 9.81 \end{matrix}] .$

微分方程式を関数ファイルとして記述し、MATLAB^® パスに保存します。

function f = cannonfodder(t,x)

f = [x(3);x(4);x(3);x(4)]; % Initial, gets f(1) and f(2) correct
nrm = norm(x(3:4)) * .02; % Norm of the velocity times constant
f(3) = -x(3)*nrm; % Horizontal acceleration
f(4) = -x(4)*nrm - 9.81; % Vertical acceleration

壁から 30 m 離れたところから pi/3 の角度で ODE の解を視覚化します。

図を生成するコード

手順 3: パターン検索を使用して解決します。

問題は、発射体が着地する壁からの距離を最大化するために、初期位置 x0(1) と初期角度 x0(2) を見つけることです。これは最大化問題なので、距離の負の値を最小化します (最大化と最小化を参照)。

patternsearch を使用してこの問題を解決するには、目的、制約、初期推定値、およびオプションを指定する必要があります。

これら 2 つのファイルは、目的関数と制約関数です。それらを MATLAB パス上のフォルダーにコピーします。

function f = cannonobjective(x)
x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];

sol = ode45(@cannonfodder,[0,15],x0);

% Find the time t when y_2(t) = 0
zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(end)]);
% Then find the x-position at that time
f = deval(sol,zerofnd,1);

f = -f; % take negative of distance for maximization



function [c,ceq] = cannonconstraint(x)

ceq = [];
x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];

sol = ode45(@cannonfodder,[0,15],x0);

if sol.y(1,end) <= 0 % Projectile never reaches wall
    c = 20 - sol.y(2,end);
else
    % Find when the projectile crosses x = 0
    zerofnd = fzero(@(r)deval(sol,r,1),[sol.x(2),sol.x(end)]);
    % Then find the height there, and subtract from 20
    c = 20 - deval(sol,zerofnd,2);
end

目的関数と制約関数が入力変数 x0 を ODE ソルバーの 4 次元初期点に設定することに注意してください。発射体が壁に当たっても、ODE ソルバーは停止しません。代わりに、制約関数は単に正になり、実行不可能な初期値を示します。

初期位置 x0(1) は 0 より上にはならず、-200 より下にするのは無駄です。（空気抵抗がなければ、最長の軌道は角度 pi/4 で –20 から始まるため、 –20 付近になるはずです。）同様に、初期角度 x0(2) は 0 未満になることはできず、 pi/2 を超えることもできません。これらの初期値から少し離れた境界を設定します。

lb = [-200;0.05];
ub = [-1;pi/2-.05];
x0 = [-30,pi/3]; % Initial guess

UseCompletePoll オプションを true に設定。これにより、より高品質なソリューションが得られ、並列計算にはこの設定が必要であるため、並列処理との直接比較が可能になります。

opts = optimoptions('patternsearch','UseCompletePoll',true);

patternsearch を呼び出してこの問題を解きます。

tic % Time the solution
[xsolution,distance,eflag,outpt] = patternsearch(@cannonobjective,x0,...
    [],[],[],[],lb,ub,@cannonconstraint,opts)
toc

Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.

xsolution =

  -28.8123    0.6095


distance =

 -125.9880


eflag =

     1


outpt = 

  struct with fields:

         function: @cannonobjective
      problemtype: 'nonlinearconstr'
       pollmethod: 'gpspositivebasis2n'
    maxconstraint: 0
     searchmethod: []
       iterations: 5
        funccount: 269
         meshsize: 8.9125e-07
         rngstate: [1×1 struct]
          message: 'Optimization terminated: mesh size less than options.MeshTolerance↵ and constraint violation is less than options.ConstraintTolerance.'

Elapsed time is 0.792152 seconds.

発射体を壁から約 29 メートル離れた場所から 0.6095 ラジアンの角度で発射すると、最も遠い距離である約 126 メートルになります。目的関数は壁までの距離の負数であるため、報告される距離は負になります。

解決策を視覚化します。

x0 = [xsolution(1);0;300*cos(xsolution(2));300*sin(xsolution(2))];

sol = ode45(@cannonfodder,[0,15],x0);
% Find the time when the projectile lands
zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(end)]);
t = linspace(0,zerofnd); % equal times for plot
xs = deval(sol,t,1); % Interpolated x values
ys = deval(sol,t,2); % Interpolated y values
plot(xs,ys)
hold on
plot([0,0],[0,20],'k') % Draw the wall
xlabel('Horizontal distance')
ylabel('Trajectory height')
legend('Trajectory','Wall','Location','NW')
ylim([0 70])
hold off

手順 4: コストのかかるサブルーチンを 2 回呼び出すことは避けてください。

目的関数と非線形制約関数は両方とも、ODE ソルバーを呼び出して値を計算します。ソルバーを 2 回呼び出すことを避けるには、同じ関数における目的と非線形制約のテクニックを使用します。runcannon ファイルはこの手法を実装しています。このファイルを MATLAB パス上のフォルダーにコピーします。

function [x,f,eflag,outpt] = runcannon(x0,opts)

if nargin == 1 % No options supplied
    opts = [];
end

xLast = []; % Last place ode solver was called
sol = []; % ODE solution structure

fun = @objfun; % The objective function, nested below
cfun = @constr; % The constraint function, nested below

lb = [-200;0.05];
ub = [-1;pi/2-.05];

% Call patternsearch
[x,f,eflag,outpt] = patternsearch(fun,x0,[],[],[],[],lb,ub,cfun,opts);

    function y = objfun(x)
        if ~isequal(x,xLast) % Check if computation is necessary
            x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];
            sol = ode45(@cannonfodder,[0,15],x0);
            xLast = x;
        end
        % Now compute objective function
        % First find when the projectile hits the ground
        zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(end)]);
        % Then compute the x-position at that time
        y = deval(sol,zerofnd,1);
        y = -y; % take negative of distance
    end

    function [c,ceq] = constr(x)
       ceq = [];
        if ~isequal(x,xLast) % Check if computation is necessary
            x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];
            sol = ode45(@cannonfodder,[0,15],x0);
            xLast = x;
        end
        % Now compute constraint functions
        % First find when the projectile crosses x = 0
        zerofnd = fzero(@(r)deval(sol,r,1),[sol.x(1),sol.x(end)]);
        % Then find the height there, and subtract from 20
        c = 20 - deval(sol,zerofnd,2);
    end

end

問題を再初期化し、runcannon の呼び出し時間を計測します。

x0 = [-30;pi/3];
tic
[xsolution,distance,eflag,outpt] = runcannon(x0,opts);
toc

Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.
Elapsed time is 0.670715 seconds.

ソルバーは以前よりも速く実行されました。ソリューションを調べると、出力は同一であることがわかります。

手順 5: 並列で計算します。

並列計算によってさらに時間を節約するようにしてください。まず、並列ワーカープールを開きます。

parpool

Starting parpool using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

ans = 

 ProcessPool with properties: 

            Connected: true
           NumWorkers: 6
              Cluster: local
        AttachedFiles: {}
    AutoAddClientPath: true
          IdleTimeout: 30 minutes (30 minutes remaining)
          SpmdEnabled: true

並列計算を使用するようにオプションを設定し、ソルバーを再実行します。

opts = optimoptions('patternsearch',opts,'UseParallel',true);
x0 = [-30;pi/3];
tic
[xsolution,distance,eflag,outpt] = runcannon(x0,opts);
toc

Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.
Elapsed time is 1.894515 seconds.

この場合、並列コンピューティングは遅くなりました。ソリューションを調べると、出力は同一であることがわかります。

手順 6.遺伝的アルゴリズムと比較してください。

遺伝的アルゴリズムを使用して問題を解決することもできます。ただし、遺伝的アルゴリズムは通常、速度が遅く、信頼性が低くなります。

同じ関数における目的と非線形制約で説明したように、ga を使用する場合、ステップ 4 で patternsearch が使用するネストされた関数手法によって時間を節約することはできません。代わりに、適切なオプションを設定して ga を並列に呼び出します。

options = optimoptions('ga','UseParallel',true);
rng default % For reproducibility
tic % Time the solution
[xsolution,distance,eflag,outpt] = ga(@cannonobjective,2,...
    [],[],[],[],lb,ub,@cannonconstraint,options)
toc

Optimization terminated: average change in the fitness value less than options.FunctionTolerance
 and constraint violation is less than options.ConstraintTolerance.

xsolution =

  -37.6217    0.4926


distance =

 -122.2181


eflag =

     1


outpt = 

  struct with fields:

      problemtype: 'nonlinearconstr'
         rngstate: [1×1 struct]
      generations: 4
        funccount: 9874
          message: 'Optimization terminated: average change in the fitness value less than options.FunctionTolerance↵ and constraint violation is less than options.ConstraintTolerance.'
    maxconstraint: 0
       hybridflag: []

Elapsed time is 12.529131 seconds.

ga ソリューションは patternsearch ソリューションほど優れていません。122 m 対 126 m。ga はより時間がかかります: 約 12 秒対 2 秒未満。patternsearch はシリアルおよびネストでさらに時間がかかり、1 秒未満です。ga を連続して実行すると、1 回のテスト実行で約 30 秒と、さらに時間がかかります。