このページは機械翻訳を使用して翻訳されました。元の英語を参照するには、ここをクリックします。
ODE を並列に最適化する
この例では、ODE のパラメータを最適化する方法を示します。
また、ODE ソリューションが目的関数と非線形制約関数の両方を返す場合に、それらの関数を 2 回計算しないようにする方法も示します。この例では、ソルバーの実行時間とソリューションの品質の観点から、patternsearch
と ga
を比較します。
並列コンピューティングを使用するには、Parallel Computing Toolbox™ ライセンスが必要です。
手順 1: 問題を定義します。
問題は、大砲の位置と角度を変えて、できるだけ壁を越えて発射物を発射することです。大砲の砲口速度は300m/sです。壁の高さは20メートルです。大砲が壁に近すぎると、急な角度で発射しなければならなくなり、発射物が十分遠くまで飛びません。大砲が壁から離れすぎている場合も、発射物は十分な距離を飛ばないことになります。
空気抵抗により発射物の速度が遅くなります。抵抗力は速度の二乗に比例し、比例定数は 0.02 です。重力は発射体に作用し、一定の 9.81 m/s2 で下向きに加速します。したがって、軌道x(t)の運動方程式は
ここで、 です。
初期位置 x0
と初期速度 xp0
は 2 次元ベクトルです。ただし、初期の高さ x0(2)
は 0 なので、初期位置はスカラー x0(1)
のみによって決まります。また、初期速度 xp0
は大きさが 300 (銃口速度) なので、初期角度 (スカラー) のみに依存します。初期角度 th
の場合、 xp0
= 300*(cos(th),sin(th))
です。したがって、最適化問題は 2 つのスカラーのみに依存するため、2 次元の問題になります。水平距離と角度を決定変数として使用します。
手順 2: ODE モデルを定式化します。
ODE ソルバーでは、モデルを 1 次システムとして定式化する必要があります。軌道ベクトル (x1(t),x2(t)) にその時間微分 (x'1(t),x'2(t)) を加えて、4 次元軌道ベクトルを形成します。この拡大ベクトルに関して、微分方程式は次のようになる。
微分方程式を関数ファイルとして記述し、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 秒と、さらに時間がかかります。