ドキュメンテーション センター

  • 評価版
  • 製品アップデート

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

2 タンク システム:連続時間 SISO システムの C MEX ファイル モデリング

この例では、C MEX モデル ファイルに基づいて IDNLGREY モデリングを実行する方法を説明します。この例では単純なシステムを使用するため、非線形の状態空間モデリングが適しています。

2 タンク システム

ここでの目的は、実験室規模の 2 タンク システムの下部タンクの水位をモデリングすることです。図 1 に概要を示します。

図 1: 2 タンク システムの概略図

入出力データ

モデル化のジョブを、使用できる入出力データの読み込みから始めます。このデータは、以下の IDNLGREY モデル構造を使用してシミュレートし、出力にノイズを追加したものです。twotankdata.mat ファイルには、サンプリング レート (Ts) 0.2 秒を使用して生成された 3000 個の入出力サンプルがあるデータセットが 1 つあります。入力 u(t) はポンプに印加された電圧 [V] で、上部タンクへの流入を生成します。この上部タンクの底部にある小さな穴から流出し、下部タンクに流れ込みます。2 タンク システムの出力 y(t) は下部タンクの水位 [m] です。タンク データを格納する IDDATA オブジェクト z を作成します。記録とドキュメンテーションのため、チャネル名と単位も指定します。このステップはオプションです。

load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'twotankdata'));
z = iddata(y, u, 0.2, 'Name', 'Two tanks');
set(z, 'InputName', 'Pump voltage', 'InputUnit', 'V',             ...
       'OutputName', 'Lower tank water level', 'OutputUnit', 'm', ...
       'Tstart', 0, 'TimeUnit', 's');

推定に使用される入出力データがプロット ウィンドウに表示されます。

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

図 2: 2 タンク システムからの入出力データ

2 タンク システムのモデル化

次のステップは、2 タンク システムを記述するモデル構造を指定することです。このためには、x1(t) と x2(t) でそれぞれ上部タンクと下部タンクの水位を示します。各タンクについて、基礎的な物理特性 (質量バランス) から、水量の変化は流入と流出の差に依存します (i = 1、2)。

  d/dt (Ai*xi(t)) = Qini(t) - Qouti(t)

ここで、Ai [m^2] はタンク i の断面積で、Qini(t) と Qouti(t) [m^3/s] は時刻 t におけるタンク i への流入とタンク i からの流出です。

上部タンクについて、流入はポンプに印加される電圧に比例すると仮定され、Qin1(t) = k*u(t) となります。上部タンクの排出孔は小さいため、Bernoulli の法則が適用され、流出は水位の平方根に比例します。正確には、次のようになります。

  Qout1(t) = a1*sqrt(2*g*x1(t))

ここで a1 は排出孔の断面積で、g は重力定数です。下部タンクでは、流入量は上部タンクからの流出量と等しくなり (Qin2(t) = Qout1(t))、流出量は Bernoulli の法則によって与えられます。

  Qout2(t) = a2*sqrt(2*g*x2(t))

ここで a2 は排出孔の断面積です。

これらの事実をまとめると、次の状態空間構造が得られます。

  d/dt x1(t) = 1/A1*(k*u(t) - a1*sqrt(2*g*x1(t)))
  d/dt x2(t) = 1/A2*(a1*sqrt(2*g*x1(t)) - a2*sqrt(2*g*x2(t)))
        y(t) = x2(t)

2 タンク C MEX モデル ファイル

次に、これらの方程式を、A1、k、a1、g、A2、a2 の 6 つのパラメーター (または定数) で C MEX ファイルに入力します。C MEX ファイルは通常、MATLAB 言語を使用して作成したファイルよりも多少複雑ですが、C MEX モデリングは一般に、特に複雑なモデルの場合、実行速度において明白な利点があります。コード記述を支援するため、テンプレート C MEX ファイルが提供されています (以下を参照してください)。多くのアプリケーションの場合、いくつかの出力を定義して、dx と y をこのテンプレートに記述するコード行を入力するだけで十分です。IDNLGREY C MEX ファイルは、常に 2 つの出力を返すように構成してください。

  dx: the right-hand side(s) of the state-space equation(s)
   y: the right-hand side(s) of the output equation(s)

また、次のように 3+Npo(+1) 入力引数を指定してください。

   t: the current time
   x: the state vector at time t ([] for static models)
   u: the input vector at time t ([] for time-series models)
   p1, p2, ..., pNpo: the individual parameters (which can be real
      scalars, column vectors or 2-dimensional matrices); Npo is here
      the number of parameter objects, which for models with scalar
      parameters coincide with the number of parameters Np
   FileArgument: optional inputs to the model file

この 2 タンク システムでは、6 つのスカラー パラメーターがあり、C MEX モデリング ファイルへの入力引数の数は 3+Npo = 3+6 = 9 になります。オプションの FileArgument がこのアプリケーションには使用されていないため、後続の 10 番目の引数は省略できます。

C MEX モデリング ファイルの作成は通常、4 つのステップで行われます。

  1. Inclusion of C-libraries and definitions of the number of outputs.
  2. Writing the function computing the right-hand side(s) of the state
     equation(s), compute_dx.
  3. Writing the function computing the right-hand side(s) of the output
     equation(s), compute_y.
  4. Writing the main interface function, which includes basic error
     checking functionality, code for creating and handling input and
     output arguments, and calls to compute_dx and compute_y.
Let us view the C MEX source file (except for some comments) for the two
tank system and based on this discuss these four items in some more
detail.

図 3: 2 タンク システムの C MEX ソース コード

1. mex.h と math.h の 2 つの C ライブラリが通常、MEX 関連の関数と数学関数にアクセスするために含まれています。出力の数も、標準 C 定義を使用してモデリング ファイルごとに宣言されます。

  /* Include libraries. */
  #include "mex.h"
  #include "math.h"
  /* Specify the number of outputs here. */
  #define NY 1

2-3. 次に、ファイル内の状態を更新する関数 compute_dx と出力を更新する関数 compute_y を検討します。この関数は両方とも引数リストを保持し、計算する出力 (dx または y) が 1 番目の引数位置にあります。これに、状態方程式および出力方程式の右辺を計算するために必要なすべての変数とパラメーターが後続します。

これらの関数での最初のステップは、以降の方程式で使用されるモデル パラメーターを展開することです。有効な変数名 (入力引数リストで使用される変数を除く) を使用して、個々のパラメーターの物理的に意味のある名前を指定できます。

C の場合と同様に、配列の最初の要素は 0 の位置に格納されます。したがって、C での dx[0] は MATLAB の dx(1) (またはスカラーの場合は単に dx) に相当し、入力 u[0] は u (または u(1))、パラメーター A1[0] は A1 に相当します。

2 タンク モデル ファイルには、平方根の計算が含まれます。これは数学 C ライブラリ math.h を含めることで可能になります。数学ライブラリはよく使用される三角関数 (sin、cos、tan、asin、acos、atan など) 指数関数 (exp)、対数関数 (log、log10)、平方根 (sqrt)、べき乗 (pow)、絶対値の計算 (fabs) を実現します。math.h ライブラリは、math.h 関数が使用されている場合は常に含めなければなりません。それ以外の場合は、省略できます。C 数値演算ライブラリについての詳細は、「非線形グレー ボックス モデルの同定に関するチュートリアル: IDNLGREY モデル ファイルの作成」を参照してください。

4. メイン インターフェイス関数は、ほぼ常に同じ内容をもち、大半のアプリケーションについて変更は不要です。原則として、変更が必要と思われるのは、compute_dx および compute_y への呼び出しが行われる場合のみです。静的システムの場合は、compute_dx への呼び出しを省くことができます。その他の場合、状態方程式と出力方程式で参照される変数とパラメーターのみを渡すことを推奨します。たとえば、1 つの状態のみが使用される 2 タンク システムの出力方程式では、入力引数リストを次のように短縮できます。

  void compute_y(double *y, double *x)

また、次のようにして compute_y を呼び出します。

  compute_y(y, x);

compute_dx と compute_y の入力引数リストを拡張して、状態の数やパラメーターの数など、インターフェイス関数で参照されるその他の変数を含めることもできます。

モデル ソース ファイルが完成したら、コンパイルしなければなりません。これは mex コマンドを使用して MATLAB コマンド プロンプトから実行できます。「help mex」を参照してください(このステップはここでは省略します)。

モデル固有の C MEX ファイルを作成する場合、IDNLGREY C MEX テンプレート ファイルをコピーして作業を開始すると便利な場合があります。このテンプレートには、スケルトン ソース コードと、特定のアプリケーション向けにコードをカスタマイズするための詳細な手順が含まれています。テンプレート ファイルの場所は、MATLAB コマンド プロンプトで次の内容を入力すると表示されます。

  fullfile(matlabroot, 'toolbox', 'ident', 'nlident', 'IDNLGREY_MODEL_TEMPLATE.c')

IDNLGREY C MEX モデル ファイルに関する詳細は、「IDNLGREY モデル ファイルの作成」の例も参照してください。

2 タンク IDNLGREY モデル オブジェクトの作成

次のステップは、2 タンク システムを記述する IDNLGREY オブジェクトを作成することです。便宜上、入力と出力に関する記録情報も設定します (名前と単位)。

FileName      = 'twotanks_c';          % File describing the model structure.
Order         = [1 1 2];               % Model orders [ny nu nx].
Parameters    = {0.5; 0.0035; 0.019; ...
                 9.81; 0.25; 0.016};   % Initial parameters.
InitialStates = [0; 0.1];              % Initial initial states.
Ts            = 0;                     % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ...
                'Name', 'Two tanks');
set(nlgr, 'InputName', 'Pump voltage', 'InputUnit', 'V',             ...
          'OutputName', 'Lower tank water level', 'OutputUnit', 'm', ...                          ...
          'TimeUnit', 's');

状態の名前と単位とモデル パラメーターに関する情報を、SETINIT および SETPAR コマンドで追加していきます。さらに、状態 x1(t) と x2(t) の両方は負に設定できないタンク水位であるため、'Minimum' プロパティで x1(0) および x2(0) が 0 以上であることも指定します。実際、すべてのモデル パラメーターはゼロより大きくなければいけないことがわかっています。そこで、すべてのパラメーターの 'Minimum' プロパティをある小さな正数 (eps(0)) に設定します。これらの設定により、推定値に関する制約が以降の推定ステップで実行されることになります (推定されたモデルは、すべての入力された制約が適用されるモデルになります)。

nlgr = setinit(nlgr, 'Name', {'Upper tank water level' 'Lower tank water level'});
nlgr = setinit(nlgr, 'Unit', {'m' 'm'});
nlgr = setinit(nlgr, 'Minimum', {0 0});   % Positive levels!
nlgr = setpar(nlgr, 'Name', {'Upper tank area'        ...
                      'Pump constant'          ...
                      'Upper tank outlet area' ...
                      'Gravity constant'       ...
                      'Lower tank area'        ...
                      'Lower tank outlet area'});
nlgr = setpar(nlgr, 'Unit', {'m^2' 'm^3/(s*V)' 'm^2'  'm/(s^2)' 'm^2' 'm^2'});
nlgr = setpar(nlgr, 'Minimum', num2cell(eps(0)*ones(6,1)));   % All parameters > 0!

2 つのタンクの断面積 (A1 と A2) は、より正確に判定できます。したがって、これらと g を定数として取り扱い、'Fixed' フィールドがすべての 6 つのパラメーターに対してコマンド GETPAR で適切に設定されていることを確認します。これらすべてによって、モデル パラメーターのうち 3 つが推定されます。

nlgr.Parameters(1).Fixed = true;
nlgr.Parameters(4).Fixed = true;
nlgr.Parameters(5).Fixed = true;
getpar(nlgr, 'Fixed')
ans = 

    [1]
    [0]
    [0]
    [1]
    [1]
    [0]

初期 2 タンク モデルの性能

自由パラメーター k、a1、a2 を推定する前に、初期パラメーター値を使用してシステムをシミュレートします。既定の微分方程式ソルバー (適応ステップ長を調整した Runge-Kutta 45 ソルバー) を使用して、絶対的および相対的許容誤差を小さな値に設定します (それぞれ 1e-6 と 1e-5)。COMPARE コマンドは、2 つの入力引数で呼び出されると、どの初期状態が 'Fixed' に定義されていても、既定の設定ですべての初期状態を推定することに注意してください。自由な初期状態のみを推定するためには、3 番目と 4 番目の入力引数を次のように指定して COMPARE を呼び出します。compare(z, nlgr, 'init', 'm')。タンク モデルの両方の初期状態は既定の設定では 'Fixed' であるため、このコマンドでは初期状態推定は実行されません。

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

compare(z, nlgr);

図 4: 初期 2 タンク モデルの実際の出力とシミュレーションの出力の比較

シミュレーションの出力と実際の出力がプロット ウィンドウに表示され、それほど一致していないことがわかります。

パラメーター推定

一致を高めるため、次に 3 つの自由パラメーターを PEM で推定します (既定の設定では、すべての初期状態の 'Fixed' フィールドは False なので、初期状態の推定は PEM へのこの呼び出しでは実行されません)。

nlgr = pem(z, nlgr, 'Display', 'Full');

2 タンク推定モデルの性能

推定モデルの性能を調査するため、そのシミュレーションが実行されます (初期状態はここで再度推定されます)。

compare(z, nlgr);

図 5: 推定 2 タンク モデルの実際の出力とシミュレーションの出力の比較

実際の出力とシミュレートした出力の一致度は高くなっています。残りの疑問は、2 タンク システムをよりシンプルな線形モデル構造を使用して正確に記述できるかどうかということです。これに答えるため、データをある標準線形モデル構造に適合させ、COMPARE を使用してこれらのモデルがタンクのダイナミクスを適切に表現しているかを確認します。

nk = delayest(z);
arx22 = arx(z, [2 2 nk]);   % Second order linear ARX model.
arx33 = arx(z, [3 3 nk]);   % Third order linear ARX model.
arx44 = arx(z, [4 4 nk]);   % Third order linear ARX model.
oe22  = oe(z,  [2 2 nk]);   % Second order linear OE model.
oe33  = oe(z,  [3 3 nk]);   % Third order linear OE model.
oe44  = oe(z,  [4 4 nk]);   % Fourth order linear OE model.
sslin = ssest(z);           % State-space model (order determined automatically)
compare(z, nlgr, 'b', arx22, 'm-', arx33, 'm:', arx44, 'm--', ...
        oe22, 'g-', oe33, 'g:', oe44, 'g--', sslin, 'r-');

図 6: 複数の推定された 2 タンク モデルの実際の出力とシミュレーションの出力の比較

比較プロットでは、線形モデルは 2 タンク システムのすべてのダイナミクスを取得することはできないことが、明確にわかります。一方で、推定された非線形 IDNLGREY モデルは、実際の出力と非常によく一致しています。さらに、IDNLGREY モデル パラメーターも実際の出力の生成に使用されたパラメーターと一致しています。IDNLGREY オブジェクトのモデル パラメーターを格納する構造体配列からパラメーター ベクトルを返すコマンド GETPVEC を使用して、表示を行います。

disp('   True      Estimated parameter vector');
ptrue = [0.5; 0.005; 0.02; 9.81; 0.25; 0.015];
fprintf('   %1.4f    %1.4f\n', [ptrue'; getpvec(nlgr)']);
   True      Estimated parameter vector
   0.5000    0.5000
   0.0050    0.0049
   0.0200    0.0200
   9.8100    9.8100
   0.2500    0.2500
   0.0150    0.0147

PE を使用して取得された予測誤差は小さく、ランダム ノイズと非常によく似ています。

figure;
pe(z, nlgr);

図 7: IDNLGREY 2 タンク モデルの推定モデルに取得した予測誤差

入力電圧が 5 から6、7、8、9、10 V に段階的に増加した場合にどうなるかも調べてみます。この作業を行うには、5 ボルトの固定オフセットから順番に複数のステップ振幅を指定して STEP を呼び出します。stepDataOptions で作成された専用のオプションセットを使用すると、ステップ応答の構成が簡単になります。

figure('Name', [nlgr.Name ': step responses']);
t = (-20:0.1:80)';
Opt = stepDataOptions('InputOffset',5,'StepAmplitude',6);
step(nlgr, t, 'b', Opt);
hold on

Opt.StepAmplitude = 7;  step(nlgr, t, 'g', Opt);
Opt.StepAmplitude = 8;  step(nlgr, t, 'r', Opt);
Opt.StepAmplitude = 9;  step(nlgr, t, 'm', Opt);
Opt.StepAmplitude = 10; step(nlgr, t, 'k', Opt);

grid on;
legend('5 -> 6 V', '5 -> 7 V', '5 -> 8 V', '5 -> 9 V', '5 -> 10 V', ...
       'Location', 'Best');

図 8: IDNLGREY 2 タンク モデルの推定モデルに取得したステップ応答

最後に、PRESENT コマンドを使用して、推定モデルの概要情報を取得します。

present(nlgr);
                                                                                   
nlgr =                                                                             
Continuous-time nonlinear grey-box model defined by 'twotanks_c' (MEX-file):       
                                                                                   
   dx/dt = F(t, u(t), x(t), p1, ..., p6)                                           
    y(t) = H(t, u(t), x(t), p1, ..., p6) + e(t)                                    
                                                                                   
 with 1 input, 2 states, 1 output, and 3 free parameters (out of 6).               
                                                                                   
 Input:                                                                            
    u(1)  Pump voltage(t) [V]                                                      
 States:                                  initial value                            
    x(1)  Upper tank water level(t) [m]   xinit@exp1     0   (fix) in [0, Inf]     
    x(2)  Lower tank water level(t) [m]   xinit@exp1   0.1   (fix) in [0, Inf]     
 Output:                                                                           
    y(1)  Lower tank water level(t) [m]                                            
 Parameters:                            value        standard dev                  
    p1   Upper tank area [m^2]                 0.5            0   (fix) in ]0, Inf]
    p2   Pump constant [m^3/(s*V)]      0.00488443    0.0259044   (est) in ]0, Inf]
    p3   Upper tank outlet area [m^2]    0.0199722   0.00645303   (est) in ]0, Inf]
    p4   Gravity constant [m/(s^2)]           9.81            0   (fix) in ]0, Inf]
    p5   Lower tank area [m^2]                0.25            0   (fix) in ]0, Inf]
    p6   Lower tank outlet area [m^2]    0.0146504    0.0776102   (est) in ]0, Inf]
                                                                                   
The model was estimated from the data set 'Two tanks', which                       
contains 3000 data samples.                                                        
Loss function 2.41439e-05 and Akaike's FPE 2.41922e-05                             
Created:       09-Aug-2013 03:03:25                                                
Last modified: 09-Aug-2013 03:03:32                                                

まとめ

この例では、次の内容について説明しました。

  1. how to use C MEX-files for IDNLGREY modeling, and
  2. provided a rather simple example where nonlinear state-space
     modeling shows good potential

その他の情報

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

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