線形グレーボックスモデルの推定

線形グレーボックスモデル構造の指定

単出力および多出力の時間領域データまたは時系列データ (出力のみ) を使用して任意の常微分方程式または差分方程式について線形の離散時間および連続時間のグレーボックスモデルを推定できます。

状態空間形式でシステム方程式を表す必要があります。"状態空間モデル" は、状態変数 "x(t)" を使用して、1 つ以上の "n" 階微分方程式ではなく、一連の 1 階微分方程式としてシステムを記述します。

グレーボックスモデリングの最初の手順では、ユーザー定義のパラメーターおよびモデルに関する情報の関数として状態空間行列を返す関数を記述します。

次の形式を使用して、ファイルで線形グレーボックスモデルを実装します。

[A,B,C,D] = myfunc(par1,par2,...,parN,Ts,aux1,aux2,...)

ここで、出力引数は状態空間行列であり、myfunc はファイルの名前です。par1,par2,...,parN はモデルの N 個のパラメーターです。各エントリはスカラー、ベクトル、または行列にすることができます。Ts はサンプル時間です。aux1,aux2,... はパラメーターおよびサンプル時間に加えて状態空間行列を計算するために myfunc で使用されるオプションの入力引数です。aux にはシステム内の補助変数が含まれます。補助変数を使用することで、ファイルを編集することなく、関数への入力でシステムパラメーターを変更します。

myfunc の内容を記述して、連続時間または離散時間の状態空間モデル、あるいはその両方をパラメーター化できます。myfunc を使用して線形グレーボックスモデルを作成すると、myfunc の出力引数の性質を指定できます。連続時間状態空間モデルの形式は以下のとおりです。

連続時間では、状態空間記述は次の形式になります。

$\begin{array}{l} \dot{x} (t) = A x (t) + B u (t) + K e (t) \\ y (t) = C x (t) + D u (t) + e (t) \\ x (0) = x 0 \end{array}$

ここで、A、B、C、および D は、パラメーター par1,par2,...,parN によってパラメーター化された行列です。ノイズ行列 K および初期状態ベクトル x0 は myfunc によってパラメーター化されません。一部のアプリケーションでは、A、B、C、および D の各行列と同様に、選択したパラメーターによってパラメーター化された数量として K および x0 を表すことができます。そのようなケースを処理するには、ODE ファイル myfunc を記述して、K および x0 を追加の出力引数として返します。

[A,B,C,D,K,x0] = myfunc(par1,par2,...,parN,Ts,aux1,aux2,...)

これにより、K および x0 は、A、B、C、および D の各行列と同様に扱われます。これらはすべて、パラメーター par1,par2,...,parN の関数です。推定時に初期状態、x0、および外乱成分 K の処理を構成するには、greyestOptions オプションセットを使用します。

離散時間では、状態空間記述は次のような似た形式になります。

$\begin{array}{l} x (k + 1) = A x (k) + B u (k) + K e (t) \\ y (k) = C x (k) + D u (k) + e (t) \\ x (0) = x 0 \end{array}$

ここで、A、B、C、および D は、今度はパラメーター par1,par2,...,parN でパラメーター化された離散時間行列です。K および x0 は直接パラメーター化されていませんが、対応する推定オプションを構成することで必要に応じて推定できます。

モデル構造をもつ関数または MEX ファイルを作成した後に、idgrey モデルオブジェクトを定義する必要があります。

グレーボックスモデルを表す関数の作成

スクリプトを開く

この例では、以下の連続時間モデルの構造を表す方法を示します。

$\begin{array}{l}
\dot x(t) = \left[ {\begin{array}{*{20}{c}}
0&1\\
0&{{\theta _1}}
\end{array}} \right]x(t) + \left[ {\begin{array}{*{20}{c}}
0\\
{{\theta _2}}
\end{array}} \right]u(t)\\
y(t) = \left[ {\begin{array}{*{20}{c}}
1&0\\
0&1
\end{array}} \right]x(t) + e(t)\\
x(0) = \left[ {\begin{array}{*{20}{c}}
{{\theta _3}}\\
0
\end{array}} \right]
\end{array}$

この方程式は電気モーターを表します。ここで、 ${y_1}(t) = {x_1}(t)$ はモーターシャフトの角度位置で、 ${y_2}(t) = {x_2}(t)$ は角速度です。パラメーター $- {\theta _1}$ はモーターの逆時定数、 $- {{\theta _2}}/{{\theta _1}}$ は入力から角速度への静的ゲインです。

モーターは t = 0 で静止していますが、その角度位置 ${\theta _3}$ は不明です。不明なパラメーターの近似ノミナル値が ${\theta _1} = - 1$ 、 ${\theta _2} = 0.25$ 、および ${\theta _3} = 0$ であるとします。この例の詳細については、System Identification: Theory for the User, Second Edition, by Lennart Ljung, Prentice Hall PTR, 1999 の状態空間モデルに関するセクションを参照してください。

連続時間状態空間モデル構造は、以下の方程式で定義されます。

$\begin{array}{l}
\dot x(t) = Fx(t) + Gu(t) + \tilde Kw(t)\\
y(t) = Hx(t) + Du(t) + w(t)\\
x(0) = x0
\end{array}$

構造化状態空間表現を使用して同じモデルを推定する場合は、Estimate Structured Continuous-Time State-Space Modelsを参照してください。

推定用にこのモデルを準備するには、以下のようにします。

次のファイルを作成して、この例のモデル構造を表します。

function [A,B,C,D,K,x0] = myfunc(par,T)
A = [0 1; 0 par(1)]; 
B = [0;par(2)];
C = eye(2);
D = zeros(2,1);
K = zeros(2,2);
x0 = [par(3);0];

MATLAB® 検索パスに含まれるようにファイルを保存します。

次の構文を使用して、myfunc ファイルに基づいてidgreyモデルオブジェクトを定義します。

par = [-1; 0.25; 0];
aux = {};
T = 0;
m = idgrey('myfunc',par,'c',aux,T);

ここで、par はすべてのユーザー定義パラメーターのベクトルを表し、ノミナル (初期) 値を含みます。この例では、すべてのスカラー値のパラメーターは par ベクトルでグループ化されます。スカラー値のパラメーターは、ODE 関数 myfunc への独立入力引数として扱うこともできました。'c' は、基となるパラメーター化が連続時間であることを指定します。aux はオプション引数を表します。myfunc にはオプション引数がないため、aux = {} を使用します。T はサンプル時間を指定します。T = 0 は連続時間モデルを示します。

推定データを読み込みます。

load dcmotordata
data = iddata(y,u,0.1);

greyestを使用して、グレーボックスパラメーター値を推定します。

m_est = greyest(data,m);

ここで、data は推定データ、m は推定初期化 idgrey モデルです。m_est は推定 idgrey モデルです。

参考

idgrey | greyest