信号伝送システム: オプションの入力引数を使用して C MEX ファイルモデリング

スクリプトを開く

この例では、IDNLGREY モデルにオプションの入力引数を設定する方法を説明します。ここでは C MEX タイプのモデルファイルに対して実行する方法を中心に説明しますが、MATLAB ファイルモデリングと類似している部分についても少し取り上げます。

ここでは、次の図で示す信号伝送システム (電信) をもとに説明します。

図 1: 信号伝送システムの概略図

信号伝送モデリング

電信の距離 x (入力電圧が印加された場所から算出) で、時刻 t で電流 i(t, x) が流れます。対応する電圧は u(t, x) で、電流と電圧の関係は 2 つの組み合わせた偏微分方程式 (PDE) で記述できます。

$- \frac{\partial u(t, x)}{\partial x} = L \frac{\partial i(t, x)}{\partial t}$

$- \frac{\partial i(t, x)}{\partial x} = C \frac{\partial u(t, x)}{\partial t}$

上記の方程式は「電信方程式」(その変種) と呼ばれることが多く、L と C はそれぞれ単位長あたりのインダクタンスとキャパシタンスです。

電信方程式は、通常の微分方程式系によって近似できますが、離散距離 0、h、2h、...、Nh での電流と電圧を考慮する場合に限ります。ここで、h は離散距離で、N は距離の合計数です。この近似を行った後、電信は互いにチェーン状に接続された、構造が同一の区間で構成されると考えられます。文献では、このタイプの近似を一般に集約と呼んでいます。

距離 x = kh (k = 0, 1, ..., N)、時刻 t での電圧と電流をそれぞれ、u_k(t) と i_k(t) とします。d/dx u(t, x) と d/dx i(t, x) を、以下の単純な差分近似で近似します。

  d u(t, x)   u_{k+1}(t) - u_k(t)
  --------- ~ -------------------           Forward approximation.
     dx               h

  d i(t, x)   i_k(t) - i_{k-1}(t)
  --------- ~ -------------------           Backward approximation.
     dx               h

x = kh

d/dx u(t, x) に前進近似を使用する理由は、i_N(t) = 0 (裸電線) として、残りの N 個の離散電流を次の N 個の微分方程式でモデル化できるようにするためです。

   d i_k(t)     1
   -------- = - -- (u_{k+1}(t) - u_k(t)) for k = 0, 1, 2, ..., N-1
      dt        Lh

同様に、u_0(t) は既知の入力信号で記述には微分方程式は不要なので、点 h、2h、...、N での d/dx i(t, x) をモデル化するためには後退近似法を使用するのが便利です。

   d u_k(t)     1
   -------- = - -- (i_k(t) - i_{k-1}(t)) for k = 1, 2, ..., N-1
      dt        Ch

   d u_N(t)     1                         1
   -------- = - -- (i_N(t) - i_{N-1}(t) = -- i_{N-1}(t)
      dt        Ch                        Ch

これによって、図 1 に示すモデル近似に達し、方程式は相互に連結されたコイルとコンデンサの数で表現されます。

次に 2*N 個の状態、x1(t) = i_0(t)、x2(t) = u_1(t)、x3(t) = i_1(t)、x4(t) = u_2(t)、...、x2N-1(t) = i_{N-1}(t)、x2N(t) = u_N(t) を導入します。また、入力 u(t) = u_0(t) 、出力を電線端部の電圧 y(t) = x2N(t) = u_N(t) とします。自明な代入により、次の状態空間モデル構造が得られます。

     x1(t) = -1/(Lh)*(x2(t) - u(t))
     x2(t) = -1/(Ch)*(x3(t) - x1(t))
     x3(t) = -1/(Lh)*(x4(t) - x2(t))
     x4(t) = -1/(Ch)*(x5(t) - x3(t))
          ...
  x2N-1(t) = -1/(Lh)*(x2N(t) - x2N-2(t))
    x2N(t) =  1/(Ch)*x2N-1(t)

      y(t) = x2N(t)

一般に、上記の数学的操作によって標準の線形状態空間モデル構造が得られ、IDNLGREY の線形版である IDGREY でうまく処理できます。ここでは IDGREY モデリングは行いませんが、オプションの入力引数を IDNLGREY で使用してモデリングの柔軟性を高める方法を説明します。

IDNLGREY 信号伝送モデルオブジェクト

柔軟性を得るため、任意の長さ L の電線をすぐに取り扱える IDNLGREY モデルオブジェクトを用意することが望まれます。モデリングの品質を高めるためには、集約ブロック N の数を変化させ、良好なシステム近似を得られるようにすると便利です。これらの要件は IDNLGREY オブジェクトの FileArgument プロパティの N と L を渡すことで処理できます。FileArgument プロパティは cell 配列でなければなりませんが、この配列には任意の種類のデータを格納できます。このアプリケーションでは、構造に N と L を指定子、長さ 1000 m の電線に対して 3 つの異なる値の N、10、30、100 を試します。次の 3 つの FileArgument がその後、IDNLGREY モデリングを実行するときに使用されます。

FileArgument10  = {struct('N', 10, 'L', 1000)};   % N = 10  --> 20 states.
FileArgument30  = {struct('N', 30, 'L', 1000)};   % N = 30  --> 60 states.
FileArgument100 = {struct('N', 100, 'L', 1000)};  % N = 100 --> 200 states.

FileArgument に含まれるデータの解析とチェックは、IDNLGREY モデルファイルで実行しなければならず、モデルファイル設計者がこの機能の実装を担当します。エラーチェックのない関数で h と N を得るには、次のコマンドを使用できます。

  N = varargin{1}{1}.N;
  h = varargin{1}{1}.L/N;

ここで FileArgument はモデルファイルに渡される最後の引数である varargin{1} に対応します。ファイル signaltransmission_m.m は信号伝送モデルを実装します。"type signaltransmission_m.m" と入力して、ファイル全体を表示します。

C MEX モデルファイルを使用するとき (後に実行) よりも状況は多少複雑です。この場合、状態 Nx の数は前もってわかりませんが、メインインターフェイス関数で計算され、入力引数として compute_dx と compute_y に渡すことができます。これらの関数の宣言は、次のようになります。

void compute_dx(double *dx, int nx, double *x, double *u, double **p,
             const mxArray *auxvar)
void compute_y(double *y, int nx, double *x)

ここで、t 変数 (方程式では使用されていません) を省き、正数 nx を両方の関数への 2 番目の引数として使用しています。さらに、出力の計算には使用されないため、compute_y の標準の末尾の 3 つの引数を削除しました。これらの変更によって、compute_dx と compute_y はメインインターフェイス関数から次のようにして呼び出されます。

  /* Call function for state derivative update. */
  compute_dx(dx, nx, x, u, p, auxvar);

  /* Call function for output update. */
  compute_y(y, nx, x);

compute_dx の最初の部分を以下に示します ("type signaltransmission_c.c" を入力すると、MATLAB® コマンドウィンドウにファイル全体が表示されます)。

  /* State equations. */
  void compute_dx(double *dx, int nx, double *x, double *u, double **p,
                  const mxArray *auxvar)
  {
     /* Declaration of model parameters and intermediate variables. */
     double *L, *C;      /* Model parameters.                  */
     double h, Lh, Ch;   /* Intermediate variables/parameters. */
     int    j;           /* Equation counter.                  */

     /* Retrieve model parameters. */
     L = p[0];   /* Inductance per unit length.  */
     C = p[1];   /* Capacitance per unit length. */

     /* Get and check FileArgument (auxvar). */
     if (mxGetNumberOfElements(auxvar) < 1) {
         mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidFileArgument",
                           "FileArgument should at least hold one element.");
     } else if (mxIsStruct(mxGetCell(auxvar, 0)) == false) {
         mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidFileArgument",
                           "FileArgument should contain a structure.");
     } else if (   (mxGetFieldNumber(mxGetCell(auxvar, 0), "N") < 0)
                || (mxGetFieldNumber(mxGetCell(auxvar, 0), "L") < 0)) {
         mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidFileArgument",
                           "FileArgument should contain a structure with fields 'N' and 'L'.");
     } else {
         /* Skip further error checking to obtain execution speed. */
         h = *mxGetPr(mxGetFieldByNumber(mxGetCell(auxvar, 0), 0,
                      mxGetFieldNumber(mxGetCell(auxvar, 0), "L"))) / (0.5*((double) nx));
     }
     Lh = -1.0/(L[0]*h);
     Ch = -1.0/(C[0]*h);
     ...

FileArgument は変数 auxvar の compute_dx に渡されることに注目してください。モデルパラメーターの宣言と取得の後 (および、h などの中間変数の宣言の後)、いくつかの外部インターフェイスルーチン (mx-routine と呼ばれる) によって auxvar の整合性がチェックされます。これらの MATLAB ルーチンによって、mxArray 変数を作成、アクセス、操作、破壊できます。この詳細については、外部インターフェイスについての MATLAB ドキュメンテーションを参照してください。最後の else 句は、すべてのチェックが合格した場合に実行されます。この場合、auxvar フィールドの値 N を使用して h の値を判定します。この h とモデルパラメーター L と C が使用され、必要なパラメーター量 Lh = -1/(L*h) と Ch = -1/(C*h) が計算されます。FileArgument についての詳細は、「非線形グレーボックスモデルの同定に関するチュートリアル: IDNLGREY モデルファイルの作成」を参照してください。

Lh と Ch を計算すると、compute_dx の 2 番目の部分は次のようになります。特に、nx が compute_dx の for loop で使用され、現在のモデルの状態数を定義する方法に注意してください。

      ...
      /* x[0]   : Current i_0(t).    */
      /* x[1]   : Voltage u_1(t).    */
      /* x[2]   : Current i_1(t).    */
      /* x[3]   : Voltage u_1(t).    */
      /* ...                         */
      /* x[Nx-2]: Current i_Nx-1(t). */
      /* x[Nx-1]: Voltage u_Nx(t).   */
      for (j = 0; j < nx; j = j+2) {
          if (j == 0) {
              /* First transmitter section. */
              dx[j]   = Lh*(x[j+1]-u[0]);
              dx[j+1] = Ch*(x[j+2]-x[j]);
          } else if (j < nx-3) {
              /* Intermediate transmitter sections. */
              dx[j]   = Lh*(x[j+1]-x[j-1]);
              dx[j+1] = Ch*(x[j+2]-x[j]);
          } else {
              /* Last transmitter section. */
              dx[j]   = Lh*(x[j+1]-x[j-1]);
              dx[j+1] = -Ch*x[j];
          }
      }
  }

出力更新関数 compute_dy は、状態更新関数よりも単純です。

  /* Output equation. */
  void compute_y(double *y, int nx, double *x)
  {
      /* y[0]: Voltage at the end of the transmitter. */
      y[0] = x[nx-1];
  }

上記の情報を 3 つの異なる IDNLGREY オブジェクト (N = 10、N = 30、N = 100 にそれぞれ 1 つずつ) に入力できる準備ができました。これらのファイルを作成する際の違いは、順序、初期状態ベクトル、使用されるオプションの入力引数だけです。

FileName        = 'signaltransmission_c';            % File describing the model structure.
Parameters      = struct('Name',    {'Inductance per unit length'    ... % Initial parameters.
                                     'Capacitance per unit length'}, ...
                         'Unit',    {'H/m' 'F/m'},                   ...
                         'Value',   {0.99e-3 0.99e-3},               ...
                         'Minimum', {eps(0) eps(0)},                 ... % L, C > 0!
                         'Maximum', {Inf Inf},                       ...
                         'Fixed',   {false false});

% A. Signal transmission model with N = 10;
Order10         = [1 1 2*FileArgument10{1}.N];       % Model orders [ny nu nx].
InitialStates10 = zeros(2*FileArgument10{1}.N, 1);   % Initial intitial states.
nlgr10 = idnlgrey(FileName, Order10, Parameters, InitialStates10, 0, ...
                'FileArgument', FileArgument10,                      ...
                'Name', '10 blocks', 'TimeUnit', 's');

% B. Signal transmission model with N = 30;
Order30         = [1 1 2*FileArgument30{1}.N];       % Model orders [ny nu nx].
InitialStates30 = zeros(2*FileArgument30{1}.N, 1);   % Initial value of intitial states.
nlgr30 = idnlgrey(FileName, Order30, Parameters, InitialStates30, 0, ...
                'FileArgument', FileArgument30,                      ...
                'Name', '30 blocks', 'TimeUnit', 's');

% C. Signal transmission model with N = 100;
Order100         = [1 1 2*FileArgument100{1}.N];     % Model orders [ny nu nx].
InitialStates100 = zeros(2*FileArgument100{1}.N, 1); % Initial value of intitial states.
nlgr100 = idnlgrey(FileName, Order100, Parameters, InitialStates100, 0, ...
                'FileArgument', FileArgument100,                        ...
                'Name', '100 blocks', 'TimeUnit', 's');

3 つの信号伝送モデルの名前と単位も、モデルのサイズが指定通り確認されると、設定されます。入力は電線に印加される電圧、出力は電線端部の電圧です。

set(nlgr10, 'InputName', 'Vin', 'InputUnit', 'V', ...
            'OutputName', 'Vout', 'OutputUnit', 'V');
set(nlgr30, 'InputName', nlgr10.InputName, 'InputUnit', nlgr10.InputUnit, ...
            'OutputName', nlgr10.OutputName, 'OutputUnit', nlgr10.OutputUnit);
set(nlgr100, 'InputName', nlgr10.InputName, 'InputUnit', nlgr10.InputUnit, ...
             'OutputName', nlgr10.OutputName, 'OutputUnit', nlgr10.OutputUnit);
nlgr10
nlgr30
nlgr100

nlgr10 =

Continuous-time nonlinear grey-box model defined by 'signaltransmission_c' (MEX-file):

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

 with 1 input(s), 20 state(s), 1 output(s), and 2 free parameter(s) (out of 2).

Name: 10 blocks

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

nlgr30 =

Continuous-time nonlinear grey-box model defined by 'signaltransmission_c' (MEX-file):

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

 with 1 input(s), 60 state(s), 1 output(s), and 2 free parameter(s) (out of 2).

Name: 30 blocks

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

nlgr100 =

Continuous-time nonlinear grey-box model defined by 'signaltransmission_c' (MEX-file):

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

 with 1 input(s), 200 state(s), 1 output(s), and 2 free parameter(s) (out of 2).

Name: 100 blocks

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

入出力データ

長さ 1000 m の信号伝送線からシミュレーションを実行した入出力データを使用できます。このデータは、上記の集約モデル構造を使用してシミュレーションを実行したものですが、より大きな N (1500) を使用しています。シミュレーションは 20 秒間行われ、サンプリングレート 0.1 秒を使用しました。使用されたモデルパラメーターは L = C = 1e-3 で、開始点はゼロ電圧線です (すべての初期状態がゼロです)。両方のモデルパラメーターは通常の信号伝送線よりも明らかに高く、このタイプのシステムで指定される転送された遅延を適切に表現するために選択されています。このデータを読み込んで、IDDATA オブジェクト z に格納します。

load('signaltransmissiondata');
z = iddata(vout, vin, 0.1, 'Name', 'Signal transmission',                ...
           'InputName', 'Vin', 'InputUnit', 'V', ...
           'OutputName', 'Vout',               ...
           'OutputUnit', 'V', 'Tstart', 0, 'TimeUnit', 's');

入出力データのプロットによって、印加された電圧から電線端部の出力電圧への伝達遅延が明確に示されます。最初の入力電圧パルス (約 1.4 秒) が約 1 秒後に出力に現れます。

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

図 2: 信号伝送システムからの入出力データ

初期信号伝送モデルの性能

これらの 3 つの初期モデルの性能はどの程度でしょうか。COMPARE を使用したモデルシミュレーションを実施して、この点を調べてみます。実行の点からは、初期状態が推定されていないため、初期状態ベクトルにゼロを選択する点が重要です。この理由は、特に nlgr100 の場合は状態の数が非常に高く (= 200)、初期状態ベクトルの推定の計算時間は非常に長くなるということです。

compare(z, nlgr100, nlgr30, nlgr10, compareOptions('InitialCondition', 'zero'));

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

一致の点では、予期されるように 3 つのモデルには大きな違いがあり、最も複雑なモデルはその他の 2 つのモデルよりも性能が高くなっています。モデル化能力の違いは、各モデルの予測誤差を確認することでよくわかります。

peOpt = peOptions('InitialCondition','zero');
e = {pe(z, nlgr100, peOpt) pe(z, nlgr30, peOpt) pe(z, nlgr10, peOpt)};
figtitle = {'nlgr100 : e@' 'nlgr30 : e@' 'nlgr10 : e@'};
figure('Name', [z.Name ': prediction errors']);
for i = 1:3
    subplot(3, 1, i);
    plot(e{i}.SamplingInstants, e{i}.OutputData, 'r');
    title(['Initial ' figtitle{i} z.OutputName{1}]);
    axis('tight');
end

図 4: 3 つの初期 IDNLGREY 信号伝送モデルで得られた予測誤差

パラメーター推定

実際のシステムの正しいパラメーター値から開始していますが、これらの値から最適なモデル一致を作成できるとは限りません。3 つの信号伝送モデルの 2 つのモデルパラメーターを NLGREYEST を使用して推定し、この点を調べてみます。これらの計算には時間がかかります。

opt = nlgreyestOptions('Display', 'on', 'EstimateCovariance', false);

nlgr100 = nlgreyest(z, nlgr100, opt);
nlgr30  = nlgreyest(z, nlgr30, opt);
nlgr10  = nlgreyest(z, nlgr10, opt);

推定された信号伝送モデルの性能

推定モデルに適用される最終予測誤差 (FPE) 基準から、最も複雑なモデルが優れていることがわかります。

fpe(nlgr100, nlgr30, nlgr10)

    0.5228    0.6771    4.7749

推定されたパラメーター値は、開始値と比較してそれほど変更されていません。以下からわかるように、あまり使用されていない 2 つのモデルのパラメーター値は実際には実パラメーター値から乖離し、最も複雑なモデルとは反対であることが興味深い点です。このようになる理由は複数あります。集約が粗い、最小化手順によって局所的最小値に対応するパラメーター値になる、などです。

disp('   True        nlgr100     nlgr30      nlgr10');
ptrue = [1e-3; 1e-3];
fprintf('   %1.7f   %1.7f   %1.7f   %1.7f\n', ...
       [ptrue'; getpvec(nlgr100)'; getpvec(nlgr30)'; getpvec(nlgr10)']);

   True        nlgr100     nlgr30      nlgr10
   0.0010000   0.0009952   0.0009759   0.0009879
   0.0010000   0.0009952   0.0009759   0.0009879

次に、COMPARE をもう一度利用して、3 つの推定された信号伝送モデルのシミュレーションを実行します。

figure
compare(z, nlgr100, nlgr30, nlgr10, compareOptions('InitialCondition', 'zero'));

図 5: 3 つの推定された信号伝送モデルの実際の出力とシミュレーションの出力の比較

3 つの例すべてで、一致は多少改善されています。パラメーター推定の前後で取得した予測誤差を比較して、小さいながらも改善を認識できます。

peOpt = peOptions('InitialCondition','zero');
e = {pe(z, nlgr100, peOpt) pe(z, nlgr30, peOpt) pe(z, nlgr10, peOpt)};
figtitle = {'nlgr100 : e@' 'nlgr30 : e@' 'nlgr10 : e@'};
figure('Name', [z.Name ': prediction errors']);
for i = 1:3
    subplot(3, 1, i);
    plot(e{i}.SamplingInstants, e{i}.OutputData, 'r');
    title(['Estimated ' figtitle{i} z.OutputName{1}]);
    axis('tight');
end

図 6: 3 つの推定された IDNLGREY 信号伝送モデルで得られた予測誤差

このチュートリアルのまとめとして、推定された信号伝送モデルの単位ステップ応答を調べます。これからわかるように、これらのすべてのモデルで伝達遅延を表現できるように思われます。ただし、これらの線形モデルの極はすべて虚数軸上、つまり安定境界上にあります。このようなモデルのシミュレーションの実行には注意が必要で、この図に示すような種類の振動ステップ応答が発生することがあることが、よく知られています。

ここでの主な問題は、この種類の電信方程式は損失をモデル化しないことです。このような効果は "抵抗" (コイルと直列) と "コンダクタンスの抵抗" (コンデンサと並列) を図 1 の各 Lh-Ch サブブロックに追加することで付加できます。これによって、モデルはより安定し、シミュレーションは簡単になります。

figure('Name', [z.Name ': step responses']);
step(nlgr100, nlgr30, nlgr10);
grid on;
legend('nlgr100', 'nlgr30', 'nlgr10', 'Location', 'NorthWest');

図 7: 3 つの推定された IDNLGREY 信号伝送モデルのステップ応答

まとめ

このチュートリアルでは、IDNLGREY モデリングを実行する際のオプションの入力引数の使用方法を最優先で説明しました。C MEX ファイルタイプのモデル化を中心として、MATLAB ファイルでの実行方法も説明しました。特定の種類のシステムについて、通常よりも柔軟なモデルファイル (使用する状態の数) を設計すると便利なことがわかります (モデルの再利用)。このチュートリアルで扱った信号伝送システムは、このカテゴリでの一例です。

信号伝送システム: オプションの入力引数を使用して C MEX ファイル モデリング