グレーボックス推定に 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)] となります。モーターの物理パラメーターに対して tau と k がどのように関連しているかについては、Ljung (1999) を参照してください。

入出力データについて

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

load('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)'});

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

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

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

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

size(nlgr)

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

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

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

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

ans = 

  struct with fields:

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


ans = 

  struct with fields:

       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 =

  2x1 cell array

    {[1]}
    {[1]}


ans =

  2x1 cell array

    {[-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(s), 2 state(s), 2 output(s), and 2 free parameter(s) (out of 2).

Name: DC-motor

Status:                                                         
Created by direct construction or transformation. Not estimated.

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

get(nlgr)

             FileName: 'dcmotor_m'
                Order: [1x1 struct]
           Parameters: [2x1 struct]
        InitialStates: [2x1 struct]
         FileArgument: {}
    SimulationOptions: [1x1 struct]
         TimeVariable: 't'
        NoiseVariance: [2x2 double]
            InputName: {'Voltage'}
            InputUnit: {'V'}
           InputGroup: [1x1 struct]
           OutputName: {2x1 cell}
           OutputUnit: {2x1 cell}
          OutputGroup: [1x1 struct]
                Notes: [0x1 string]
             UserData: []
                 Name: 'DC-motor'
                   Ts: 0
             TimeUnit: 'seconds'
               Report: [1x1 idresults.nlgreyest]

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

パラメーター tau および k を推定する前に、既定の微分方程式ソルバー (適応ステップ長を調整したルンゲ・クッタ 45 ソルバー) を使用してシステムの出力を推定したパラメーターでシミュレートします。シミュレーションオプションは、"SimulationOptions" モデルプロパティを使って指定されます。

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

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

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

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

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

compare(z, nlgr);

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

パラメーター推定

パラメーターと初期状態を、非線形グレーボックスモデルの予測誤差最小化法である nlgreyest を使用して推定します。推定の進捗状況表示を選択するなどの推定オプションは、オプションセット "nlgreyestOptions" を使用して指定します。

nlgr = setinit(nlgr, 'Fixed', {false false}); % Estimate the initial states.
opt = nlgreyestOptions('Display', 'on');
nlgr = nlgreyest(z, nlgr, opt);

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

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

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

nlgr.Report
fprintf('\n\nThe search termination condition:\n')
nlgr.Report.Termination

ans = 

         Status: 'Estimated using NLGREYEST'
         Method: 'Solver: ode45; Search: lsqnonlin'
            Fit: [1x1 struct]
     Parameters: [1x1 struct]
    OptionsUsed: [1x1 idoptions.nlgreyest]
      RandState: []
       DataUsed: [1x1 struct]
    Termination: [1x1 struct]



The search termination condition:

ans = 

  struct with fields:

                 WhyStop: 'Change in cost was less than the specified tolerance.'
              Iterations: 5
    FirstOrderOptimality: 1.4014e-04
                FcnCount: 6
               Algorithm: 'trust-region-reflective'

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.残差 (「レフトオーバー」) を確認します。

残差は、モデルで説明されないままの部分を示し、小さければモデルの品質は良好です。残差間の相関を表示するには resid コマンドを使用します。プロットの最初の列は、2 つの出力の残差における自己相関を示します。2 列目は、これらの残差と入力 "Voltage" との相互相関を示します。相関は許容範囲 (青い領域) 内にあります。

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 番目のモデルパラメーターがシステムダイナミクスの説明に重要であることを示しています。小さなノイズ分散 (複数出力システムの共分散) 要素も、モデルが推定データを正しく取得できることを示す良い指標です。

getcov(nlgr)
nlgr.NoiseVariance

ans =

   1.0e-04 *

    0.1573    0.0021
    0.0021    0.0008


ans =

    0.0010   -0.0000
   -0.0000    0.0110

推定されたモデルの詳細については、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(s), 2 state(s), 2 output(s), and 2 free parameter(s) (out of 2).

 Inputs:
    u(1)  Voltage(t) [V]
 States:                              Initial value
    x(1)  Angular position(t) [rad]     xinit@exp1   0.0302675   (estimated) in [-Inf, Inf]
    x(2)  Angular velocity(t) [rad/s]   xinit@exp1   -0.133777   (estimated) in [-Inf, Inf]
 Outputs:
    y(1)  Angular position(t) [rad]
    y(2)  Angular velocity(t) [rad/s]
 Parameters:                        Value Standard Deviation
    p1   Time-constant [s]           0.243649       0.00396671   (estimated) in [-Inf, Inf]
    p2   Static gain [rad/(V*s)]     0.249644      0.000284486   (estimated) in [-Inf, Inf]

Name: DC-motor

Status:                                                                         
Termination condition: Change in cost was less than the specified tolerance..   
Number of iterations: 5, Number of function evaluations: 6                      
                                                                                
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "DC-motor".
Fit to estimation data: [98.34;84.47]%                                          
FPE: 0.001096, MSE: 0.1187                                                      
More information in model's "Report" property.

まとめ

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

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

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

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

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

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

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