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

IDNLGREY モデル ファイルの作成

この例では、非線形グレー ボックス モデルの ODE ファイルを MATLAB ファイルおよび C MEX ファイルとして作成する方法を説明します。

グレー ボックス モデリングは、より包括的なモデル化手順を伴う点でブラックボックス モデリングとは概念的に異なります。IDNLGREY (非線形グレーボックス モデル オブジェクト、IDGREY の非線形に相当) の場合、この手順では ODE ファイル (モデル ファイル) を作成します。ODE ファイルでは、物理的な第一原理モデリングを通じて到達する状態方程式と出力方程式の右辺を指定します。この例では、このファイルを MATLAB ファイルまたは C MEX ファイルとして実装する場合の一般的な方法を説明します。

IDNLGREY モデル ファイル

IDNLGREY は、以下の明示的な状態空間形式 (いわゆる出力誤差 (OE) 形式 - ノイズ e(t) がモデル構造の出力に加法的にのみ影響することからこのように名付けられました) で記述された非線形モデル構造におけるパラメーターと初期状態の推定をサポートします。

xn(t) = F(t, x(t), u(t), p1, ..., pNpo); x(0) = X0;

y(t) = H(t, x(t), u(t), p1, ..., pNpo) + e(t)

離散時間構造の場合は xn(t) = x(T+Ts) (Ts はサンプル時間)、連続時間構造の場合は xn(t) = d/dt x(t) です。さらに、F(.) と H(.) は任意の線形関数または非線形関数で、それぞれ Nx (状態の数) 個または Ny (出力の数) 個の成分をもちます。モデル パラメーター (p1, ..., pNpo のうちどれでも) と初期状態ベクトル X(0) を推定できます。以下の点に注目してください。

  1. 時系列モデリング、つまり外力入力信号 u(t) のないモデリング、および

  2. 静的モデリング、つまり状態 x(t) のないモデリング

は、IDNLGREY がサポートする 2 つの特殊なケースです。(この 2 つのモデリング カテゴリの例は、チュートリアル idnlgreydemo3 および idnlgreydemo5 を参照してください)。

IDNLGREY のモデリングでは最初に必ず、状態の更新と出力の計算の方法を指定する MATLAB または C MEX のモデル ファイルを実装します。さらに重要なこととして、以下の入力引数および出力引数 (この形式は MATLAB タイプと C MEX タイプのどちらのモデル ファイルでも必要です) を使って定義したモデル ファイル MODFILENAME.m または MODFILENAME.c を作成しなければなりません。

[dx, y] = MODFILENAME(t, x, u, p1, p2, ..., pNpo, FileArgument)

MODFILENAME には MATLAB または C MEX ファイルの任意のユーザー定義名 (twotanks_m.m や pendulum_c.c のような) を指定できます。このファイルは 2 つの出力を返すように定義しなければなりません。

  • dx: 状態空間方程式の右辺 (Nx 実数エントリをもつ列ベクトル、静的モデルの場合は [])

  • y: 出力方程式の右辺 (Ny 実数エントリをもつ列ベクトル)

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

  • t: 現在の時刻

  • x: 時間 t における状態ベクトル (静的モデルの場合は [])

  • u: 時間 t における入力ベクトル (時系列モデルの場合は [])

  • p1, p2, ..., pNpo: 個々のパラメーター (実数スカラー、列ベクトルまたは 2 次元行列)。Npo は、ここではパラメーター オブジェクトの数で、スカラー パラメーターをもつモデルの場合はパラメーターの数 Np と一致

  • FileArgument: モデル ファイルに対するオプション入力

以降、MATLAB 言語または C-MEX ファイルを使ったモデルの作成を中心に説明を進めていきます。ただし、IDNLGREY は P ファイル (MATLAB コマンド "pcode" を使って取得される保護された MATLAB ファイル) と関数ハンドルもサポートしています。実際、C MEX モデル ファイルだけでなく Fortran MEX ファイルも使用できます。後者の詳細は、MATLAB の『External Interfaces』ドキュメンテーションを参照してください。

どのような種類のモデル ファイルを実装するべきでしょうか。この疑問への答えはモデルの用途によって異なります。

MATLAB 言語を使った実装 (作成されるのは *.m ファイル)には、いくつか顕著な利点があります。第 1 に、時間のかかる低水準のプログラミングを回避し、より集中してモデル化を行うことができます。第 2 に、MATLAB とそのツールボックスに用意されている関数は、どれもモデル ファイルの中で直接使用できます。第 3 に、この種のファイルはサイズが小さく、何も変更しなくても組み込みの MATLAB エラー チェックがすべて自動的に実行されます。さらに、これはコードのコンパイルなしで取得できます。

C MEX モデリングはより複雑で、C プログラミング言語の基本的な知識が必要です。C MEX モデル ファイルの主な利点は実行速度の向上です。一般的なアドバイスとして、作成したモデルを何度も使用する場合、大きなデータセットを使用する場合、および/またはモデル構造に多くの計算が含まれる場合には、C MEX モデリングを追求することをお勧めします。多くの場合、最初は MATLAB ファイルを使用して、後で対応する C MEX ファイルに切り替えるという方法が有効です。

MATLAB 言語を使って作成された IDNLGREY モデル ファイル

上記の理由で、次は MATLAB ファイルのモデル化に進み、例として 2 タンク システムを説明する非線形 2 次モデル構造を使用します。モデル化の詳細は idnlgreydemo2 を参照してください。twotanks_m.m の内容は次のとおりです。

type twotanks_m.m
function [dx, y] = twotanks_m(t, x, u, A1, k, a1, g, A2, a2, varargin)
%TWOTANKS_M  A two tank system.

%   Copyright 2005-2006 The MathWorks, Inc.

% Output equation.
y = x(2);                                              % Water level, lower tank.

% State equations.
dx = [1/A1*(k*u(1)-a1*sqrt(2*g*x(1)));             ... % Water level, upper tank.
      1/A2*(a1*sqrt(2*g*x(1))-a2*sqrt(2*g*x(2)))   ... % Water level, lower tank.
     ];

関数ヘッダーには、必須の入力引数 t、x、および u があり、その後に 6 個のスカラー モデル パラメーター A1、k、a1、g、A2、および a2 が続きます。MATLAB ファイルの場合、オプションのモデル ファイル入力引数 FileArgument をサポートするために、最後の入力引数は常に varargin でなければなりません。IDNLGREY モデル オブジェクトでは、FileArgument はあらゆる種類のデータを保持できる cell 配列として保存されます。ここでは FileArgument の最初の要素は varargin{1}{1} を通じてアクセスされます。

変数とパラメーターは MATLAB の標準的な方法で参照されます。最初の状態は x(1) で 2 番目の状態は x(2)、入力は u(1)(スカラーの場合は u のみ)、スカラー パラメーターは各自の名前 (A1、k、a1、g、A2、および a2) を通じてアクセスされます。ベクトル パラメーターと行列パラメーターの個々の要素は、それぞれ P(i) (P という名前のベクトル パラメーターの要素 i) および P(i, j) (P という名前の行列パラメーターの行 i 列 j にある要素) としてアクセスされます。

IDNLGREY C MEX モデル ファイル

C MEX モデル ファイルの作成は MATLAB モデル ファイルの作成より複雑です。この手順を簡素化するために、利用可能な IDNLGREY C MEX モデル テンプレートを MODFILENAME.c にコピーすることをお勧めします。このテンプレートには、スケルトン ソース コードと、特定のアプリケーション向けにコードをカスタマイズするための詳細な手順が含まれています。テンプレート ファイルの場所は、MATLAB のコマンド プロンプトで次のように入力するとわかります。

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

2 つのタンクの例では、このテンプレートが twotanks_c.c にコピーされました。最初のいくつかの変更と設定の後 (以下で説明します)、状態方程式と出力方程式が入力された結果、C MEX ソース コードは以下のようになりました。

type twotanks_c.c
/*   Copyright 2005-2015 The MathWorks, Inc. */
/*   Written by Peter Lindskog. */

/* Include libraries. */
#include "mex.h"
#include <math.h>

/* Specify the number of outputs here. */
#define NY 1

/* State equations. */
void compute_dx(double *dx, double t, double *x, double *u, double **p,
                const mxArray *auxvar)
{
    /* Retrieve model parameters. */
    double *A1, *k, *a1, *g, *A2, *a2;
    A1 = p[0];   /* Upper tank area.        */
    k  = p[1];   /* Pump constant.          */
    a1 = p[2];   /* Upper tank outlet area. */
    g  = p[3];   /* Gravity constant.       */
    A2 = p[4];   /* Lower tank area.        */
    a2 = p[5];   /* Lower tank outlet area. */
    
    /* x[0]: Water level, upper tank. */
    /* x[1]: Water level, lower tank. */
    dx[0] = 1/A1[0]*(k[0]*u[0]-a1[0]*sqrt(2*g[0]*x[0]));
    dx[1] = 1/A2[0]*(a1[0]*sqrt(2*g[0]*x[0])-a2[0]*sqrt(2*g[0]*x[1]));
}

/* Output equation. */
void compute_y(double *y, double t, double *x, double *u, double **p,
               const mxArray *auxvar)
{
    /* y[0]: Water level, lower tank. */
    y[0] = x[1];
}



/*----------------------------------------------------------------------- *
   DO NOT MODIFY THE CODE BELOW UNLESS YOU NEED TO PASS ADDITIONAL
   INFORMATION TO COMPUTE_DX AND COMPUTE_Y
 
   To add extra arguments to compute_dx and compute_y (e.g., size
   information), modify the definitions above and calls below.
 *-----------------------------------------------------------------------*/

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
    /* Declaration of input and output arguments. */
    double *x, *u, **p, *dx, *y, *t;
    int     i, np;
    size_t  nu, nx;
    const mxArray *auxvar = NULL; /* Cell array of additional data. */
    
    if (nrhs < 3) {
        mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidSyntax",
        "At least 3 inputs expected (t, u, x).");
    }
    
    /* Determine if auxiliary variables were passed as last input.  */
    if ((nrhs > 3) && (mxIsCell(prhs[nrhs-1]))) {
        /* Auxiliary variables were passed as input. */
        auxvar = prhs[nrhs-1];
        np = nrhs - 4; /* Number of parameters (could be 0). */
    } else {
        /* Auxiliary variables were not passed. */
        np = nrhs - 3; /* Number of parameters. */
    }
    
    /* Determine number of inputs and states. */
    nx = mxGetNumberOfElements(prhs[1]); /* Number of states. */
    nu = mxGetNumberOfElements(prhs[2]); /* Number of inputs. */
    
    /* Obtain double data pointers from mxArrays. */
    t = mxGetPr(prhs[0]);  /* Current time value (scalar). */
    x = mxGetPr(prhs[1]);  /* States at time t. */
    u = mxGetPr(prhs[2]);  /* Inputs at time t. */
    
    p = mxCalloc(np, sizeof(double*));
    for (i = 0; i < np; i++) {
        p[i] = mxGetPr(prhs[3+i]); /* Parameter arrays. */
    }
    
    /* Create matrix for the return arguments. */
    plhs[0] = mxCreateDoubleMatrix(nx, 1, mxREAL);
    plhs[1] = mxCreateDoubleMatrix(NY, 1, mxREAL);
    dx      = mxGetPr(plhs[0]); /* State derivative values. */
    y       = mxGetPr(plhs[1]); /* Output values. */
    
    /*
      Call the state and output update functions.
      
      Note: You may also pass other inputs that you might need,
      such as number of states (nx) and number of parameters (np).
      You may also omit unused inputs (such as auxvar).
      
      For example, you may want to use orders nx and nu, but not time (t)
      or auxiliary data (auxvar). You may write these functions as:
          compute_dx(dx, nx, nu, x, u, p);
          compute_y(y, nx, nu, x, u, p);
    */
    
    /* Call function for state derivative update. */
    compute_dx(dx, t[0], x, u, p, auxvar);
    
    /* Call function for output update. */
    compute_y(y, t[0], x, u, p, auxvar);
    
    /* Clean up. */
    mxFree(p);
}

このファイルの内容を見ていきましょう。最初に気付くのは、C MEX モデル ファイルを作成する作業を以下の 4 つのサブステップに分けられるということです (最後のステップはオプション)。

  1. C ライブラリを含めて出力数を定義。

  2. 状態方程式の右辺を計算する関数 compute_dx を作成。

  3. 出力方程式の右辺を計算する関数 compute_y を作成。

  4. オプションで、基本エラー チェック機能、入出力引数を作成および処理するコード、compute_dx および compute_y に対する呼び出しを含むメイン インターフェイス関数を更新。

これらのサブステップを詳しく説明する前に、C プログラミング言語のいくつかの一般的な特徴について簡単に述べておきます。

  1. 高精度の変数 (IDNLGREY オブジェクトのすべての入力、状態、出力およびパラメーター) を、データ型が "double" になるように定義する必要があります。

  2. 変数名またはパラメーター名の直前に配置された単項の * 演算子は、いわゆる逆参照演算子です。C 宣言 "double *A1;" は、A1 が double 変数に対するポインターであることを指定します。ポインターの構成は、C 内の概念であり、常に簡単に理解できるとは限りません。幸い、compute_y と compute_dx の出力/入力変数の宣言に変更がなく、アンパックされたモデル パラメーターがすべて * を含めて内部で宣言される場合、IDNLGREY モデリングの観点からは、ポインターについて詳細を知る必要はありません。

  3. compute_y と compute_dx はいずれもまず宣言されて実装され、その後メイン インターフェイス関数で呼び出されます。宣言で、キーワード "void" は、値が何も返されないことを明示的に示します。

C プログラミング言語の詳細は、以下の文献を参照してください。

B.W. Kernighan and D. Ritchie. The C Programming Language, 2nd

edition, Prentice Hall, 1988.

最初のサブステップでは、まず C ライブラリ "mex.h" (必須) および "math.h" (高度な計算の場合は必須) を含めます。出力の数も、標準 C 定義を使用してモデリング ファイルごとに宣言されます。

/* Include libraries. */

#include "mex.h"

#include "math.h"

/* Specify the number of outputs here. */

#define NY 1

必要に応じて、上記以外の C ライブラリをさらに含めることもできます。

状態方程式または出力方程式に三角関数や平方根関数などの高度な計算が含まれている場合には、"math.h" ライブラリを含めなければなりません。"math.h" に含まれている関数と MATLAB でそれらに対応する関数のリストを以下に示します。

C-function MATLAB function

========================================

sin, cos, tan sin, cos, tan

asin, acos, atan asin, acos, atan

sinh, cosh, tanh sinh, cosh, tanh

exp, log, log10 exp, log, log10

pow(x, y) x^y

sqrt sqrt

fabs abs

MATLAB 関数の方が対応する C 関数より柔軟である点に注意してください。たとえば、MATLAB は複素数を扱いますが、C 関数は扱いません。

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

パラメーターはすべて配列 P に格納されます。compute_dx および compute_y の最初のステップは、その後の方程式で使用されるパラメーターをアンパックして名前を付けることです。twotanks_c.c において、compute_dx は 6 個のパラメーター変数を宣言します。それぞれの値は以下のように決定されます。

/* Retrieve model parameters. */

double *A1, *k, *a1, *g, *A2, *a2;

A1 = p[0]; /* Upper tank area. */

k = p[1]; /* Pump constant. */

a1 = p[2]; /* Upper tank outlet area. */

g = p[3]; /* Gravity constant. */

A2 = p[4]; /* Lower tank area. */

a2 = p[5]; /* Lower tank outlet area. */

それに対して compute_y は出力の計算にパラメーターを必要としないため、取得されるモデル パラメーターはありません。

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

上記の例ではスカラー パラメーターのみを使用しているので、パラメーターの総数 Np はパラメーター オブジェクトの数 Npo と等しくなります。ベクトルまたは行列パラメーターがモデルに含まれている場合は、Npo < Np となります。

スカラー パラメーターは P[0] (MATLAB ファイルでは P(1) または単に P) として参照され、ベクトル要素 i:th は P[i-1] (MATLAB ファイルでは P(i)) として参照されます。C MEX モデル ファイルに渡される行列は、列が互いの上に明確な順序で積み上がっていくという意味で異なります。したがって、P が 2 行 2 列の行列である場合、P(1, 1) は P[0] として参照され、P(2, 1) は P[1]、P(1, 2) は P[2]、P(2, 2) は P[3] として参照されます。スカラー、ベクトルおよび行列パラメーターが使用されている例については、「非線形グレー ボックス モデルの同定に関するチュートリアル: 産業用 3 自由度ロボット: ベクトル/行列パラメーターを使用した MIMO システムの C MEX ファイル モデリング」(idnlgreydemo8) を参照してください。

状態を更新する関数と出力を更新する関数も、単にパラメーターを取得して式の右辺を計算する以外の計算を含んでいる可能性があります。実行速度を上げるために、後続の式で値が何回も利用される中間変数を宣言して使用するという方法があります。これをよく示している例として、上記のロボットに関するチュートリアル (idnlgreydemo8) があります。

compute_dx と compute_y もオプションの FileArgument を処理できます。FileArgument データは変数 auxvar を使ってこれらの関数に渡されるため、FileArgument (cell 配列) の最初の成分は以下の式で取得できます。

mxArray* auxvar1 = mxGetCell(auxvar, 0);

ここで、mxArray は MATLAB 定義のデータ型で、C MEX ファイルと MATLAB の間でのデータの交換を可能にします。したがって、auxvar1 はあらゆるデータを格納できます。auxvar1 の解析、チェック、および使用は、これらの関数の中でのみ行われなければなりません。この機能を実装するかどうかはモデル ファイルの設計者しだいです。mxArrays を扱う関数の詳細は、MATLAB の『External Interfaces』ドキュメンテーションを参照してください。C MEX モデル ファイルのオプションの引数の使用例については、「非線形グレー ボックス モデルの同定に関するチュートリアル: 信号伝送システム: オプションの入力引数を使用して C MEX ファイル モデリング」を参照してください。

メイン インターフェイス関数は、ほぼ常に同じ内容をもち、大半のアプリケーションについて変更は不要です。原則として、変更が必要と思われるのは、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 の入力引数リストを拡張し、インターフェイス関数の中で推定される変数をさらに含めることもできます。次が計算され、それを渡すことができます。整数変数 nu (入力の数)、nx (状態の数)、および np (ここではパラメーター オブジェクトの数)。たとえばチュートリアル idnlgreydemo6 で研究するモデルでは、nx が compute_y に渡されます。

完成した C MEX モデル ファイルは IDNLGREY モデリングに使用する前にコンパイルしなければなりません。コンパイルは MATLAB コマンド ラインから以下のコマンドを使って簡単に実行できます。

mex MODFILENAME.c

mex コマンドははじめて使用する前に設定しなければなりません。それも MATLAB コマンド ラインから以下のコマンドを使って実行します。

mex -setup

IDNLGREY モデル オブジェクト

すぐに実行できるモデル ファイルが完成すれば、シミュレーション、パラメーターの推定などの実行対象となる IDNLGREY モデル オブジェクトを作成するのは簡単です。その例として、2 タンク システムを説明する 2 つの異なる IDNLGREY モデル オブジェクトを作成します。一方は MATLAB で作成したモデル ファイルを使用して、もう一方は上で詳述した C MEX ファイルを使用して作成します (C MEX モデル ファイルが既にコンパイル済みである点に注意してください)。

Order         = [1 1 2];               % Model orders [ny nu nx].
Parameters    = [0.5; 0.003; 0.019; ...
                 9.81; 0.25; 0.016];   % Initial parameter vector.
InitialStates = [0; 0.1];              % Initial values of initial states.
nlgr_m    = idnlgrey('twotanks_m', Order, Parameters, InitialStates, 0)
nlgr_m =
Continuous-time nonlinear grey-box model defined by 'twotanks_m' (MATLAB 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(s), 2 state(s), 1 output(s), and 6 free parameter(s) (out of 6).

Status:                                                         
Created by direct construction or transformation. Not estimated.
nlgr_cmex = idnlgrey('twotanks_c', Order, Parameters, InitialStates, 0)
nlgr_cmex =
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(s), 2 state(s), 1 output(s), and 6 free parameter(s) (out of 6).

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

まとめ

このチュートリアルでは MATLAB および C MEX を用いて IDNLGREY モデル ファイルを作成する方法を説明してきました。最後に、現在用意されている IDNLGREY モデル ファイルとそれらが使用されているチュートリアル/ケース スタディのリストを示して、このプレゼンテーションのまとめとします。さらなる比較を簡単にするために、MATLAB (FILENAME_m.m という形式のファイル名) モデル ファイルと C MEX モデル ファイル (FILENAME_c.c という形式のファイル名) の両方を一覧にし、「チュートリアル/ケース スタディ」の列で、どちらのタイプのモデリング アプローチがそのチュートリアルまたはケース スタディで使用されているかを示しました。

Tutorial/Case study MATLAB file C MEX-file

======================================================================

idnlgreydemo1 (MATLAB) dcmotor_m.m dcmotor_c.c

idnlgreydemo2 (C MEX) twotanks_m.m twotanks_c.c

idnlgreydemo3 (MATLAB) preys_m.m preys_c.c

(C MEX) predprey1_m.m predprey1_c.c

(C MEX) predprey2_m.m predprey2_c.c

idnlgreydemo4 (MATLAB) narendrali_m.m narendrali_c.c

idnlgreydemo5 (MATLAB) friction_m.m friction_c.c

idnlgreydemo6 (C MEX) signaltransmission_m.m signaltransmission_c.c

idnlgreydemo7 (C MEX) twobodies_m.m twobodies_c.c

idnlgreydemo8 (C MEX) robot_m.m robot_c.c

idnlgreydemo9 (MATLAB) cstr_m.m cstr_c.c

idnlgreydemo10 (MATLAB) pendulum_m.m pendulum_c.c

idnlgreydemo11 (C MEX) vehicle_m.m vehicle_c.c

idnlgreydemo12 (C MEX) aero_m.m aero_c.c

idnlgreydemo13 (C MEX) robotarm_m.m robotarm_c.c

これらのモデル ファイルの内容は、MATLAB コマンド ウィンドウで コマンド "type FILENAME_m.m" または "type FILENAME_c.c" を使って表示できます。モデル ファイルはすべて以下の MATLAB コマンドが返すディレクトリに収められています。

fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'examples')