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

並列で常微分方程式を最適化する

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

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

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

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

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

• 手順 3: patternsearch を使用して解決します。

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

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

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

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

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

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

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

ここで $$v(t)\text{}=dx(t)/dt.$$

初期位置 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}{dt}x(t)=\left[\begin{array}{c}{x}_3(t)\\ x_4(t)\\ -.02\Vert \left(x_3(t),x_4(t)\right)\Vert x_3(t)\\ -.02\Vert \left(x_3(t),x_4(t)\right)\Vert x_4(t)-9.81\end{array}\right].$$

微分方程式を関数ファイルとして記述し、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 の解を視覚化します。

 図を生成するコード

x0 = [-30;0;300*cos(pi/3);300*sin(pi/3)];
sol = ode45(@cannonfodder,[0,10],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 120])
hold off

手順 3: patternsearch を使用して解決します。

問題は、発射物が壁から着地する距離を最大化するために、初期位置 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 秒と、さらに時間がかかります。

参考

トピック

• 同じ関数における目的と非線形制約
• 並列計算
• 非線形制約によるサロゲート最適化