標準的振子: アルゴリズム関係の問題

ライブスクリプトを開く

この例では、推定アルゴリズムの選択が非線形グレーボックスモデルの推定にどのような影響を及ぼすかを説明します。ここでは、非線形振子システム (図 1) で生成されたデータを使用します。ここでは特に、微分方程式ソルバーが結果にどのような影響を及ぼすかを説明します。

図 1: 標準的振子の概略図

出力データ

出力データ (時系列データ) を読み込むことからモデル化作業を始めます。データには 1 つの出力 y があり、これは振子の角度位置 [rad] です。角度は、振子が下向きの時はゼロで、反時計回りに増加していきます。データ点のサンプル (シミュレーションを実行したもの) は 1001 あり、サンプル時間は 0.1 秒です。振子は一定の重力の影響を受けますが、振子の動作はその他の外力の影響は受けません。この状況を調べるため、IDDATA オブジェクトを作成します。

load pendulumdata
z = iddata(y, [], 0.1, 'Name', 'Pendulum');
z.OutputName = 'Pendulum position';
z.OutputUnit = 'rad';
z.Tstart = 0;
z.TimeUnit = 's';

振子の角度位置 (出力) をプロットウィンドウに表示します。

figure('Name', [z.Name ': output data']);
plot(z);

図 2: 振子の角度位置 (出力)

振子のモデル化

次のステップは、振子システムのモデル構造を指定することです。このダイナミクスは多くの書籍や資料で取り上げられ、理解されています。振子の角度位置 [rad] として x1(t)、角速度 [rad/s] として x2(t) を選択すると、次の種類の状態空間構造を簡単に設定できます。

  d/dt x1(t) = x2(t)
  d/dt x2(t) = -(g/l)*sin(x1(t)) - (b/(m*l^2))*x2(t)

        y(t) = x1(t)

パラメーター (または定数) を次のように指定します。

  g - the gravity constant [m/s^2]
  l - the length of the rod of the pendulum [m]
  b - viscous friction coefficient [Nms/rad]
  m - the mass of the bob of the pendulum [kg]

この情報を、次の内容の MATLAB ファイル pendulum_m.m に入力します。

  function [dx, y] = pendulum_m(t, x, u, g, l, b, m, varargin)
  %PENDULUM_M  A pendulum system.

  % Output equation.
  y = x(1);                                   % Angular position.

  % State equations.
  dx = [x(2);                             ... % Angular position.
        -(g/l)*sin(x(1))-b/(m*l^2)*x(2)   ... % Angular velocity.
       ];

次のステップは、振子を記述する汎用 IDNLGREY オブジェクトを作成することです。また、入力、状態、出力、パラメーターの名前と単位に関する情報も入力します。物理的な現実性を確保するため、すべてのモデルパラメーターは正でなければなりません。これはすべてのパラメーターの "Minimum" プロパティを、MATLAB® で認識できる最小の正の値である eps(0) に設定することで実現されます。

FileName      = 'pendulum_m';        % File describing the model structure.
Order         = [1 0 2];             % Model orders [ny nu nx].
Parameters    = [9.81; 1; 0.2; 1];   % Initial parameters.
InitialStates = [1; 0];              % Initial initial states.
Ts            = 0;                   % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts);
nlgr.OutputName = 'Pendulum position';
nlgr.OutputUnit = 'rad';
nlgr.TimeUnit = 's';
nlgr = setinit(nlgr, 'Name', {'Pendulum position' 'Pendulum velocity'});
nlgr = setinit(nlgr, 'Unit', {'rad' 'rad/s'});
nlgr = setpar(nlgr, 'Name', {'Gravity constant' 'Length' ...
                      'Friction coefficient' 'Mass'});
nlgr = setpar(nlgr, 'Unit', {'m/s^2' 'm' 'Nms/rad' 'kg'});
nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0) eps(0)});   % All parameters > 0.

初期の 3 つの振子モデルの性能

パラメーターの推定を始める前に、予測したパラメーター値でシステムをシミュレートします。使用できる 3 つのソルバー、固定ステップ長の前進オイラー (ode1)、適応ステップ長の Runge-Kutta 23 (ode23)、適応ステップ長の Runge-Kutta 45 (ode45) について実行します。これらの 3 つのソルバーを使用して得られた出力が、プロットウィンドウに表示されます。

% A. Model computed with first-order Euler forward ODE solver.
nlgref = nlgr;
nlgref.SimulationOptions.Solver = 'ode1';        % Euler forward.
nlgref.SimulationOptions.FixedStep = z.Ts*0.1;   % Step size.

% B. Model computed with adaptive Runge-Kutta 23 ODE solver.
nlgrrk23 = nlgr;
nlgrrk23.SimulationOptions.Solver = 'ode23';     % Runge-Kutta 23.

% C. Model computed with adaptive Runge-Kutta 45 ODE solver.
nlgrrk45 = nlgr;
nlgrrk45.SimulationOptions.Solver = 'ode45';     % Runge-Kutta 45.

compare(z, nlgref, nlgrrk23, nlgrrk45, 1, ...
   compareOptions('InitialCondition', 'model'));

図 3: 3 つの初期振子モデルの実際の出力とシミュレーションの出力の比較

見てわかるように、前進オイラー法 (既定の設定で使用されるよりもグリッドは細かい) の結果は、Runge-Kutta ソルバーからの結果と大幅に異なります。この場合、前進オイラーソルバーからは良好な結果 (モデル一致の面から) が得られ、Runge-Kutta ソルバーからの出力は多少制限されています。ただし、これは事実とは多少異なっており、後で説明しますが、粘性摩擦係数の開始値 b は実際の 2 倍になっています。

パラメーター推定

重力定数 g、棒の長さ l、振子玉の質量 m は容易に測定でき、または推定しなくてもテーブルから取得できます。しかし、摩擦係数にはこれはあてはまらず、通常は推定しなければなりません。したがって、以前の 3 つのパラメーターを g = 9.81、l = 1、m = 1と固定して、b のみを自由パラメーターとします。

nlgref = setpar(nlgref, 'Fixed', {true true false true});
nlgrrk23 = setpar(nlgrrk23, 'Fixed', {true true false true});
nlgrrk45 = setpar(nlgrrk45, 'Fixed', {true true false true});

次に、3 つの微分方程式ソルバーを使用して b を推定します。前進オイラー法に基づくモデル構造から始めます。

opt = nlgreyestOptions('Display', 'on');
tef = clock;
nlgref = nlgreyest(z, nlgref, opt);   % Perform parameter estimation.
tef = etime(clock, tef);

次に、Runge-Kutta 23 ソルバー (ode23) に基づくモデルに進みます。

trk23 = clock;
nlgrrk23 = nlgreyest(z, nlgrrk23, opt);   % Perform parameter estimation.
trk23 = etime(clock, trk23);

最後に、Runge-Kutta 45 ソルバー (ode45) を使用します。

trk45 = clock;
nlgrrk45 = nlgreyest(z, nlgrrk45, opt);   % Perform parameter estimation.
trk45 = etime(clock, trk45);

推定される 3 つの振子モデルの性能

これらの 3 つのソルバーを使用した結果を、以下にまとめます。b の実際の値は 0.1 で、ode 45 から得られます。ode23 からは 0.1 に非常に近い値が返され、ode1 は 0.1 から遠い値が返されます。

disp('           Estimation time   Estimated b value');

           Estimation time   Estimated b value

fprintf('   ode1 :  %3.1f               %1.3f\n', tef, nlgref.Parameters(3).Value);

   ode1 :  7.0               0.194

fprintf('   ode23:  %3.1f               %1.3f\n', trk23, nlgrrk23.Parameters(3).Value);

   ode23:  2.1               0.093

fprintf('   ode45:  %3.1f               %1.3f\n', trk45, nlgrrk45.Parameters(3).Value);

   ode45:  2.7               0.100

ただし、これによってシミュレーションの出力が実出力とは異なるということにはなりません。異なる微分方程式ソルバーから生じる誤差は、通常推定手順で考慮されているからです。3 つの推定された振子モデルをシミュレートすることで、この事実がすぐにわかります。

compare(z, nlgref, nlgrrk23, nlgrrk45, 1, ...
   compareOptions('InitialCondition', 'model'));

図 4: 3 つの推定された振子モデルの実際の出力とシミュレーションの出力の比較

この図から予期されたとおり、これらのモデルの最終予測誤差 (FPE) 基準値も互いに近いものです。

fpe(nlgref, nlgrrk23, nlgrrk45);

   1.0e-03 *

    0.1609    0.1078    0.1056

これに基づき、3 つのすべてのモデルで振子のダイナミクスを表現でき、前進オイラー法に基づくモデルは基礎となる物理現象をあまり反映していないことがわかります。「物理」的な性能を向上させるには、ソルバーのステップ長を減らすほかありませんが、これによって計算時間も大幅に長くなります。この実験では、精度と計算速度の両方の点で、Runge-Kutta 45 ソルバーがノンスティッフ問題については最適なソルバーであることが示されています。

まとめ

Runge-Kutta 45 ソルバー (ode45) は比較的高速で高品質の結果を返すことが多く、この理由により、IDNLGREY における既定の微分方程式ソルバーとして選択されています。「help idnlgrey.SimulationOptions」と入力すると、使用可能なソルバーについての詳細が提示され、「help nlgreyestOptions」と入力すると IDNLGREY モデリングに使用できるさまざまな推定オプションについての詳細が提示されます。