ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

グレー ボックス推定に MATLAB ファイルを使用した非線形ダイナミクスの表現

この例では、非線形グレー ボックス モデルを構成、推定、分析する方法を説明します。

非線形グレー ボックス (idnlgrey) モデルは、連続または離散時間で非線形状態空間構造によって記述されるシステムのパラメーター推定に適しています。idgrey (線形グレー ボックス モデル) と idnlgrey の両方のオブジェクトを使用して線形システムをモデル化できます。ただし、非線形ダイナミクスの表現に使用できるのは idnlgrey のみです。idgrey を使用した線形グレー ボックスのモデル化については、「System Identification Toolbox™ を使用した構造化およびユーザー定義モデルの構築」を参照してください。

モデルについて

この例では、idnlgrey オブジェクトを使用して線形 DC モーターのダイナミクスをモデル化します。

図 1: DC モーターの概略図

外乱を無視して y(1) をモーターの角度位置 [rad]、y(2) を角速度 [rad/s] として選択する場合、次の形式の線形状態空間構造を設定できます (原典については、Ljung, L 著『System Identification: Theory for the User, Upper Saddle River』NJ、Prentice-Hall PTR、1999、2nd ed.、p. 95-97 を参照してください)。

   d         | 0      1   |        |   0   |
   -- x(t) = |            | x(t) + |       | u(t)
   dt        | 0   -1/tau |        | k/tau |
             | 1   0 |
      y(t) = |       | x(t)
             | 0   1 |

tau はモーターの時定数 [s]、k は角速度への入力からの静的ゲイン [rad/(V*s)] となります。モーターの物理パラメーターに対して tauk がどのように関連しているかについては、Ljung (1999) を参照してください。

入出力データについて

1.DC モーター データを読み込みます。

load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'dcmotordata'));

2.推定データを iddata オブジェクトとして表現します。

z = iddata(y, u, 0.1, 'Name', 'DC-motor');

3.入力信号名と出力信号名、開始時間、時間単位を指定します。

z.InputName = 'Voltage';
z.InputUnit =  'V';
z.OutputName = {'Angular position', 'Angular velocity'};
z.OutputUnit = {'rad', 'rad/s'};
z.Tstart = 0;
z.TimeUnit = 's';

4.データをプロットします。

データは 2 つのプロット ウィンドウに表示されます。

figure('Name', [z.Name ': Voltage input -> Angular position output']);
plot(z(:, 1, 1));   % Plot first input-output pair (Voltage -> Angular position).
figure('Name', [z.Name ': Voltage input -> Angular velocity output']);
plot(z(:, 2, 1));   % Plot second input-output pair (Voltage -> Angular velocity).

図 2: DC モーターの入出力データ

DC モーターの線形モデリング

1.DC モーター構造を関数で表現します。

この例では、MATLAB® ファイルを使用しますが、C MEX ファイル (計算速度の向上のため)、P ファイルまたは関数ハンドルを使用することもできます。詳細は、「IDNLGREY モデル ファイルの作成」を参照してください。

DC モーター関数は dcmotor_m.m と呼ばれ、以下のようになっています。

  function [dx, y] = dcmotor_m(t, x, u, tau, k, varargin)
  % Output equations.
  y = [x(1);                         ... % Angular position.
       x(2)                          ... % Angular velocity.
      ];
  % State equations.
  dx = [x(2);                        ... % Angular velocity.
        -(1/tau)*x(2)+(k/tau)*u(1)   ... % Angular acceleration.
       ];

ファイルは常に、以下のものを返す構造にしなければなりません。

出力引数:

  • dx は、連続時間の場合は状態の微係数のベクトルで、離散時間の場合は状態更新値です。

  • y は出力方程式です。

入力引数:

  • 最初の 3 つの入力引数は、次のようにしなければなりません。t (時間)、x (状態ベクトル、静的システムの場合は [])、u (入力ベクトル、時系列の場合は [])。

  • パラメーターの順序付きリストが後に続きます。パラメーターにはスカラー、列ベクトル、2 次元行列を指定できます。

  • 補助入力引数の varargin

2.DC モーター ダイナミクスを idnlgrey オブジェクトを使用して表現します。

モデルは、状態方程式を使用して入力が出力を生成する方法を示します。

FileName      = 'dcmotor_m';       % File describing the model structure.
Order         = [2 1 2];           % Model orders [ny nu nx].
Parameters    = [1; 0.28];         % Initial parameters. Np = 2.
InitialStates = [0; 0];            % Initial initial states.
Ts            = 0;                 % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ...
                'Name', 'DC-motor');

実際には、出力に影響する外乱があります。idnlgrey モデルは外乱を明示的にはモデル化するものではありませんが、出力に追加されることを仮定しています。したがって、idnlgrey モデルは出力エラー (OE) モデルと同等です。ノイズ モデルがなければ、過去の出力が将来の出力の予測に影響することがなくなるため、予測区間 k の出力の予測はシミュレーションの出力と一致します。

3.入力と出力の名前と単位を指定します。

set(nlgr, 'InputName', 'Voltage', 'InputUnit', 'V',               ...
          'OutputName', {'Angular position', 'Angular velocity'}, ...
          'OutputUnit', {'rad', 'rad/s'},                         ...
          'TimeUnit', 's');

4.初期状態の名前と単位およびパラメーターを指定します。

nlgr = setinit(nlgr, 'Name', {'Angular position' 'Angular velocity'});
nlgr = setinit(nlgr, 'Unit', {'rad' 'rad/s'});
nlgr = setpar(nlgr, 'Name', {'Time-constant' 'Static gain'});
nlgr = setpar(nlgr, 'Unit', {'s' 'rad/(V*s)'});

setinitsetpar を使用して、値、最小値、最大値、推定ステータスをすべての初期状態またはパラメーターに同時に割り当てることもできます。

5.初期モデルを表示します。

a.モデルの基本情報を取得します。

DC モーターには 2 つの (初期) 状態と 2 つのモデル パラメーターがあります。

size(nlgr)
Nolinear grey-box model with 2 outputs, 1 inputs, 2 states and 2 parameters (2 free).

b.初期状態とパラメーターを表示します。

初期状態とパラメーターは両方とも構造体配列です。フィールドでは個々の初期状態またはパラメーターのプロパティを指定します。詳細については、idprops idnlgrey InitialStates および idprops idnlgrey Parameters を入力します。

nlgr.InitialStates(1)
nlgr.Parameters(2)
ans = 

       Name: 'Angular position'
       Unit: 'rad'
      Value: 0
    Minimum: -Inf
    Maximum: Inf
      Fixed: 1


ans = 

       Name: 'Static gain'
       Unit: 'rad/(V*s)'
      Value: 0.2800
    Minimum: -Inf
    Maximum: Inf
      Fixed: 0

c.すべての初期状態またはモデル パラメーターについての情報を 1 回の呼び出しで取得します。

たとえば、固定された (推定されていない) 初期状態とすべてのモデル パラメーターの最小値の情報を取得します。

getinit(nlgr, 'Fixed')
getpar(nlgr, 'Min')
ans = 

    [1]
    [1]


ans = 

    [-Inf]
    [-Inf]

d.オブジェクトの基本情報を取得します。

nlgr
nlgr =
Continuous-time nonlinear grey-box model defined by 'dcmotor_m' (MATLAB file):

   dx/dt = F(t, u(t), x(t), p1, p2)
    y(t) = H(t, u(t), x(t), p1, p2) + e(t)

 with 1 input, 2 states, 2 outputs, and 2 free parameters (out of 2).

get を使用して、モデル プロパティの詳細情報を取得します。idnlgrey オブジェクトはパラメトリック線形モデル オブジェクトの多数のプロパティを共有しています。

get(nlgr)
            FileName: 'dcmotor_m'
               Order: [1x1 struct]
          Parameters: [2x1 struct]
       InitialStates: [2x1 struct]
        FileArgument: {}
    CovarianceMatrix: 'estimate'
      EstimationInfo: [1x1 struct]
        TimeVariable: 't'
       NoiseVariance: [2x2 double]
           Algorithm: [1x1 struct]
                  Ts: 0
            TimeUnit: 'seconds'
           InputName: {'Voltage'}
           InputUnit: {'V'}
          InputGroup: [1x1 struct]
          OutputName: {2x1 cell}
          OutputUnit: {2x1 cell}
         OutputGroup: [1x1 struct]
                Name: 'DC-motor'
               Notes: {}
            UserData: []

初期 DC モーター モデルの性能評価

パラメーター tau および k を推定する前に、既定の微分方程式ソルバー (適応ステップ長を調整した Runge-Kutta 45 ソルバー) を使用してシステムの出力を推測したパラメーターでシミュレートします。

1.絶対許容誤差と相対許容誤差を小さな値に設定します (それぞれ 1e-61e-5)。

nlgr.Algorithm.SimulationOptions.AbsTol = 1e-6;
nlgr.Algorithm.SimulationOptions.RelTol = 1e-5;

2.シミュレートした出力と測定されたデータを比較します。

compare は 1 つまたは複数のモデルの測定された出力とシミュレーションの出力の両方を表示しますが、predict は同じ入力引数で呼び出され、シミュレーションの出力を表示します。

シミュレートした出力はプロット ウィンドウに示されます。

compare(z, nlgr);

図 3: 初期 DC モーター モデルの測定出力とシミュレーションの出力の比較

パラメーター推定

pem を使用してパラメーターと初期状態を推定します (予測誤差法)。

nlgr = setinit(nlgr, 'Fixed', {false false});   % Estimate the initial state.
nlgr = pem(z, nlgr, 'Display', 'Full');

推定 DC モーター モデルの性能評価

1.推定プロセスに関する情報を確認します。

この情報は、idnlgrey オブジェクトの EstimationInfo プロパティに保存されています。プロパティには、ソルバーや検索方法などのモデルの推定方法、データセット、推定の終了理由についての情報も含まれています。

nlgr.EstimationInfo
ans = 

             Status: 'Estimated model (PEM)'
             Method: 'Solver: ode45; Search: lsqnonlin'
            LossFcn: 0.0011
                FPE: 0.0011
           DataName: 'DC-motor'
         DataLength: 400
             DataTs: {[0.1000]}
         DataDomain: 'Time'
    DataInterSample: {'zoh'}
            WhyStop: 'Change in cost was less than the specified tolerance'
         UpdateNorm: []
    LastImprovement: []
         Iterations: 5
       InitialGuess: [1x1 struct]
            Warning: ''
     EstimationTime: 25.6100

2.シミュレーションの出力と測定された出力を比較して、モデルの品質を評価します。

一致度は 98% と 84% で、推定されたモデルは DC モーターのダイナミクスを良好に表しているといえます。

compare(z, nlgr);

図 4: 推定された IDNLGREY DC モーター モデルの測定出力とシミュレーションの出力の比較

3.idnlgrey モデルの性能を 2 次 ARX モデルと比較します。

na = [2 2; 2 2];
nb = [2; 2];
nk = [1; 1];
dcarx = arx(z, [na nb nk]);
compare(z, nlgr, dcarx);

図 5: 推定された IDNLGREY および ARX DC モーター モデルの測定出力とシミュレーションの出力の比較

4.予測誤差を確認します。

取得された予測誤差は小さく、ゼロを中心としています (バイアスなし)。

pe(z, nlgr);

図 6: IDNLGREY DC モーターの推定モデルで取得した予測誤差

5.残差 (「レフトオーバー」) を確認します。

残差は、モデルで説明されないままの部分を示し、小さければモデルの品質は良好です。以下の 2 行のコードを実行して、残差プロットを生成します。任意のキーを押して、別のプロットに進みます。

  figure('Name', [nlgr.Name ': residuals of estimated model']);
  resid(z, nlgr);

図 7: IDNLGREY DC モーターの推定モデルで取得した残差

6.ステップ応答をプロットします。

単位入力ステップで、ランプ型の動作を示す角度位置と、一定レベルで安定する角速度が生成されます。

figure('Name', [nlgr.Name ': step response of estimated model']);
step(nlgr);

図 8: IDNLGREY DC モーターの推定モデルで取得したステップ応答

7.モデルの共分散を調べます。

推定モデルの品質は、推定共分散行列と推定ノイズ分散を確認することで、ある程度評価できます。共分散行列の (i, i) 対角要素の「小さな」値は、選択したモデル構造を使用する際に、i 番目のモデル パラメーターがシステム ダイナミクスの説明に重要であることを示しています。小さなノイズ分散 (複数出力システムの共分散) 要素も、モデルが推定データを正しく取得できることを示す良い指標です。

nlgr.CovarianceMatrix
nlgr.NoiseVariance
ans =

   1.0e-04 *

    0.1521    0.0015
    0.0015    0.0007


ans =

    0.0099   -0.0004
   -0.0004    0.1094

推定されたモデルの詳細は、present を使用して初期状態および推定されたパラメーター値と、パラメーターの推定される不確かさ (標準偏差) を表示します。

present(nlgr);
                                                                                     
nlgr =                                                                               
Continuous-time nonlinear grey-box model defined by 'dcmotor_m' (MATLAB file):       
                                                                                     
   dx/dt = F(t, u(t), x(t), p1, p2)                                                  
    y(t) = H(t, u(t), x(t), p1, p2) + e(t)                                           
                                                                                     
 with 1 input, 2 states, 2 outputs, and 2 free parameters (out of 2).                
                                                                                     
 Input:                                                                              
    u(1)  Voltage(t) [V]                                                             
 States:                                initial value                                
    x(1)  Angular position(t) [rad]     xinit@exp1   0.0302675   (est) in [-Inf, Inf]
    x(2)  Angular velocity(t) [rad/s]   xinit@exp1   -0.133777   (est) in [-Inf, Inf]
 Outputs:                                                                            
    y(1)  Angular position(t) [rad]                                                  
    y(2)  Angular velocity(t) [rad/s]                                                
 Parameters:                       value      standard dev                           
    p1   Time-constant [s]         0.243649    0.00390034   (est) in [-Inf, Inf]     
    p2   Static gain [rad/(V*s)]   0.249644   0.000272168   (est) in [-Inf, Inf]     
                                                                                     
The model was estimated from the data set 'DC-motor', which                          
contains 400 data samples.                                                           
Loss function 0.00107459 and Akaike's FPE 0.00108534                                 
Created:       21-Jul-2014 19:58:48                                                  
Last modified: 21-Jul-2014 19:59:25                                                  

まとめ

この例では、非線形グレー ボックス モデル化を実行する基本ツールを説明しました。その他の非線形グレー ボックスの例を参照してください。

  • 非線形連続時間および離散時間、時系列および静的モデルの構築など、より高度なモデリング状況に非線形グレー ボックス モデルを使用する

  • C MEX モデル ファイルの作成と使用

  • 非スカラー パラメーターの処理

  • 特定のアルゴリズムの選択による影響

System Identification Toolbox を使った動的システムの同定についての詳細は、System Identification Toolbox 製品の情報ページを参照してください。

この情報は役に立ちましたか?