Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

arx

ARX、ARIX、AR、または ARI モデルのパラメーターの推定

説明

AR モデルまたは ARX モデルの推定

sys = arx(tt,[na nb nk]) は、timetable tt の変数に含まれるデータを使用して、ARX モデルまたは AR idpoly モデル sys のパラメーターを推定します。最初の Nu 個の変数が入力として使用され、次の Ny 個の変数が出力として使用されます。Nu と Ny はそれぞれ nbna の次元から決定されます。

AR モデルの場合は、入力信号がないため、sys = arx(tt,na) を使用します。この場合は、最初の Ny 個の変数を使用してモデルが適合されます。

arx は、最小二乗法と、[na nb nk] で指定された多項式の次数を使用して推定を行います。モデルのプロパティには、共分散 (パラメーターの不確かさ) および推定データと測定データの適合度が含まれます。

tt から特定の入力チャネルと出力チャネルを選択するには、名前と値の構文を使用して 'InputName''OutputName' を対応する timetable 変数名に設定します。

sys = arx(u,y,[na nb nk]) は、コンマ区切りの行列 u,y 内の時間領域入力信号と出力信号を使用します。データのサンプル時間は 1 秒と仮定されます。サンプル時間を変更するには、名前と値の構文を使用して Ts を設定します。

sys = arx(data,[na nb nk]) は、データ オブジェクト data 内の時間領域データまたは周波数領域データを使用します。この構文は特に、周波数領域データまたは周波数応答データを使用してモデルを推定する場合、またはデータ オブジェクトが提供する追加情報 (データのサンプル時間、実験のラベル付けなど) を利用する場合に使用します。

追加のオプションの指定

sys = arx(___,Name,Value) は、1 つ以上の名前と値のペアの引数を使用して追加のオプションを指定します。たとえば、名前と値のペアの引数 'IntegrateNoise',1 は、非定常外乱があるシステムに役立つ ARIX 構造または ARI 構造のモデルを推定します。この構文では、前述の任意の入力引数の組み合わせで使用できます。

sys = arx(___,opt) は、オプション セット opt を使用して推定オプションを指定します。

推定される初期条件を返す

[sys,ic] = arx(___) は、推定される初期条件を initialCondition オブジェクトとして返します。この構文は、モデルの応答を同じ推定入力データを使用してシミュレートまたは予測し、その応答を同じ推定出力データと比較する場合に使用します。初期条件を組み込むことで、シミュレーションの最初の部分における一致が高まります。

すべて折りたたむ

指定した ARX モデルに基づいて出力データを生成し、その出力データを使用してモデルを推定します。

ARX 構造をもつ多項式モデル sys0 を指定します。モデルには 1 サンプルの入力遅延が含まれており、B 多項式の先頭のゼロで表現されています。

A = [1  -1.5  0.7];
B = [0 1 0.5];
sys0 = idpoly(A,B);

ランダムなバイナリ ノイズを含む測定入力信号 u と正規分布ノイズを含む誤差信号 e を生成します。これらの信号を使用して、sys0 の測定出力信号 y をシミュレートします。

u = iddata([],idinput(300,'rbs'));
e = iddata([],randn(300,1));
y = sim(sys0,[u e]);

yu を単一の iddata オブジェクト z に結合します。z から元のモデルと同じ多項式の次数および入力遅延を使用して新しい ARX モデルを推定します。

z = [y,u];
sys = arx(z,[2 2 1])
sys =
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)
  A(z) = 1 - 1.524 z^-1 + 0.7134 z^-2              
                                                   
  B(z) = z^-1 + 0.4748 z^-2                        
                                                   
Sample time: 1 seconds
  
Parameterization:
   Polynomial orders:   na=2   nb=2   nk=1
   Number of free coefficients: 4
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using ARX on time domain data "sys0".  
Fit to estimation data: 81.36% (prediction focus)
FPE: 1.025, MSE: 0.9846                          
 

出力に、推定の他の詳細と一緒に、推定パラメーターを含む多項式が表示されます。Status の下の Fit to estimation data から、推定モデルの 1 ステップ先の予測精度が 80% を超えていることがわかります。

関数 arx を使用して時系列 AR モデルを推定します。AR モデルには測定入力はありません。

ノイズがある時系列 tt9 を含む timetable にあるデータを読み込みます。

load sdata9 tt9

[na nb nk]na の次数のみを指定して、4 次 AR モデルを推定します。

sys = arx(tt9,4);

推定された A 多項式のパラメーターとデータに対する推定の適合度を調べます。

param = sys.Report.Parameters.ParVector
param = 4×1

   -0.7923
   -0.4780
   -0.0921
    0.4698

fit = sys.Report.Fit.FitPercent
fit = 79.4835

ARIX モデルのパラメーターを推定します。ARIX モデルは、積分ノイズを含む ARX モデルです。

ARX 構造をもつ多項式モデル sys0 を指定します。モデルには 1 サンプルの入力遅延が含まれており、B の先頭のゼロで表現されています。

A = [1  -1.5  0.7];
B = [0 1 0.5];
sys0 = idpoly(A,B);

ランダムなバイナリ入力信号 u と正規分布誤差信号 e を使用して、sys0 の出力信号をシミュレートします。

u = iddata([],idinput(300,'rbs'));
e = iddata([],randn(300,1));
y = sim(sys0,[u e]);

出力信号を積分し、結果 yiiddata オブジェクト zi に格納します。

yi = iddata(cumsum(y.y),[]);
zi = [yi,u];

zi から ARIX モデルを推定します。名前と値のペアの引数 'IntegrateNoise'true に設定します。

sys = arx(zi,[2 2 1],'IntegrateNoise',true);

5 ステップの予測を使用してモデルの出力を予測し、結果を yi と比較します。

compare(zi,sys,5)

Figure contains an axes object. The axes object with ylabel y1 contains 2 objects of type line. These objects represent Validation data (y1), sys: 76.25%.

arxRegulを使用して正則化定数を自動で決定し、その値を次数 50 の FIR モデルの推定に使用します。

lambdaR の値を取得します。

load regularizationExampleData eData;
orders = [0 50 0];
[lambda,R] = arxRegul(eData,orders);

返された lambdaR の値を正則化 ARX モデルの推定に使用します。

opt = arxOptions;
opt.Regularization.Lambda = lambda;
opt.Regularization.R = R;
sys = arx(eData,orders,opt);

データを読み込みます。

load iddata1ic z1i

2 次 ARX モデル sys を推定し、ic で初期条件を返します。

na = 2;
nb = 2;
nk = 1;
[sys,ic] = arx(z1i,[na nb nk]);
ic
ic = 
  initialCondition with properties:

     A: [2x2 double]
    X0: [2x1 double]
     C: [0 2]
    Ts: 0.1000

ic は、sys の自由応答を状態空間形式で X0 の初期状態ベクトルにカプセル化する initialCondition オブジェクトです。ic は、sysz1i 入力信号でシミュレートし、その応答を z1i 出力信号と比較する場合に組み込むことができます。

入力引数

すべて折りたたむ

推定データ。一定間隔の時間ベクトルを使用する timetable として指定します。tt には入力チャネルと出力チャネルを表す変数が含まれます。複数実験データの場合、tt は長さが Ne の timetable の cell 配列です。Ne は実験数です。

指定された多項式の次数の次元から、推定に使用する入力チャネルと出力チャネルの数が判定されます。入力/出力チャネルの選択は、'InputName''OutputName' の名前と値の引数が指定されているかどうかによります。

  • 'InputName''OutputName' が指定されていない場合、tt の最初の Nu 個の変数が入力として使用され、tt の次の Ny 個の変数が出力として使用されます。

  • 'InputName''OutputName' が指定されている場合、指定された変数が使用されます。指定された入力名と出力名の数は、Nu および Ny と一致しなければなりません。

  • 入力がない時系列モデルを推定できる関数では、'InputName' を指定する必要はありません。

推定データ型の取り扱いの詳細については、Data Domains and Data Types in System Identification Toolboxを参照してください。

推定データ。SISO システムの場合は、等間隔でサンプリングされた入力および出力時間領域信号値を含む、Ns 行 1 列の実数値行列のコンマ区切りのペアとして指定します。ここで、Ns はサンプルの数です。

MIMO システムの場合は、u,y を、以下の次元をもつ入力行列と出力行列のペアとして指定します。

  • u — Ns 行 Nu 列の行列。Nu は入力の数です。

  • y — Ns 行 Ny 列の行列。Ny は出力の数です。

複数実験データの場合は、u,y を 1 行 Ne 列の cell 配列のペアとして指定します。Ne は実験数です。すべての実験のサンプル時間が一致していなければなりません。

出力のみを含み、入力を含まない時系列データの場合は、[],y を指定します。

制限

  • 行列ベースのデータは、周波数領域データからの推定をサポートしません。iddata オブジェクトや idfrd オブジェクトなどのデータ オブジェクト (data を参照) を使用する必要があります。

推定データ型の取り扱いの詳細については、Data Domains and Data Types in System Identification Toolboxを参照してください。

推定データ。iddata オブジェクト、frd (Control System Toolbox) オブジェクト、または idfrd 周波数応答オブジェクトとして指定します。AR および ARI の時系列モデルの場合、data の入力チャネルは空でなければなりません。

モデルの多項式の次数と遅延。1 行 3 列のベクトル、または行列 [na nb nk] のベクトルとして指定します。多項式の次数は、その多項式で推定する係数の数と等しくなります。

AR または ARI の時系列モデルの場合、入力はないため、[na nb nk] をスカラー na に設定します。例については、AR モデルを参照してください。

Ny 個の出力と Nu 個の入力をもつモデルの場合は次のようになります。

  • na は多項式 A(q) の次数です。Ny 行 Ny 列の非負の整数の行列として指定します。

  • nb は多項式 B(q) + 1 の次数です。Ny 行 Nu 列の非負の整数の行列として指定します。

  • nk は入出力遅延です。伝達遅延とも呼ばれます。Ny 行 Nu 列の非負の整数の行列として指定します。nk は、ARX モデルでは B 多項式の先頭に含まれる固定のゼロで表現されます。

    たとえば、伝達遅延がない状態の sys.b[5 6] であるとします。

    • sys.b + 1 は 2 次多項式であるため、nb = 2 です。

    • 伝達遅延として nk = 3 を指定します。この遅延を指定すると、sys.b の先頭にゼロが 3 つ追加され、sys.b[0 0 0 5 6] となります。nb は 2 のままです。

    • これらの係数は、多項式 B(q) = 5 q-3 + 6q-4 を表します。

    伝達遅延は、名前と値のペアの引数 'IODelay' を使用して実装することもできます。

.

例: arx(data,[2 1 1]) は、iddata オブジェクトから、1 サンプルの入力遅延を含む 1 つの入力チャネルをもつ 2 次 ARX モデルを計算します。

ARX モデルの同定の推定オプション。arOptions オプション セットとして指定します。opt で指定するオプションには以下が含まれます。

  • 初期条件の処理 — このオプションは、周波数領域データにのみ使用します。時間領域データの場合は、非測定信号が予測子で必要にならないように信号がシフトされます。

  • 入力データと出力データのオフセット — これらのオプションは、推定時に時間領域データからオフセットを削除するために使用します。

  • 正則化 — このオプションは、推定プロセスにおけるバイアスと分散の誤差のトレードオフを制御するために使用します。

詳細については、arxOptions を参照してください。例については、正則化を含む ARX モデルを参照してください。

名前と値の引数

引数のオプションのペアを Name1=Value1,...,NameN=ValueN として指定します。Name は引数名、Value は対応する値です。名前と値の引数は他の引数の後に表示されなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

例: 'IntegrateNoise',true は、ノイズ チャネルに積分器を追加します。

入力チャネル名。string、文字ベクトル、string 配列、または文字ベクトルの cell 配列として指定します。

データ ソースに timetable を使用している場合は、InputName の名前は timetable 変数のサブセットでなければなりません。

例: sys = arx(tt,__,'InputName',["u1" "u2"]) は、推定で使用する timetable tt からの入力チャネルとして変数 u1 および u2 を選択します。

出力チャネル名。string、文字ベクトル、string 配列、または文字ベクトルの cell 配列として指定します。

データ ソースに timetable を使用している場合は、OutputName の名前は timetable 変数のサブセットでなければなりません。

例: sys = arx(tt,__,'OutputName',["y1" "y3"]) は、推定で使用する timetable tt からの出力チャネルとして変数 y1 および y3 を選択します。

サンプル時間。'Ts' とサンプル時間 (TimeUnit で指定された単位) で構成されるコンマ区切りのペアとして指定します。行列ベースのデータ (u,y) を使用するときに、仮定された 1 秒のサンプル時間以外のサンプル時間が必要な場合は Ts を指定する必要があります。

timetable tt のデータのサンプル時間を取得するには、timetable プロパティ tt.Properties.Timestep を使用します。

例: arx(umat1,ymat1,___,'Ts',0.08) は、0.08 秒のサンプル時間でモデルを計算します。

入力遅延。サンプル時間の整数倍として表されます。'InputDelay' と次のいずれかで構成されるコンマ区切りのペアとして指定します。

  • Nu 行 1 列のベクトル (Nu は入力の数) — 各エントリは、対応する入力チャネルの入力遅延を表す数値です。

  • スカラー値 — すべての入力チャネルに同じ遅延を適用します。

例: arx(data,[2 1 3],'InputDelay',1) は、3 サンプルの入力遅延を含む 1 つの入力チャネルをもつ 2 次 ARX モデルを推定します。

各入出力ペアの伝達遅延。サンプル時間の整数倍として表されます。'IODelay' と次のいずれかで構成されるコンマ区切りのペアとして指定します。

  • Ny 行 Nu 列の行列 (Ny は出力の数、Nu は入力の数) — 各エントリは、対応する入出力ペアの伝達遅延を表す整数値です。

  • スカラー値 — すべての入出力ペアに同じ遅延を適用します。この方法は、入出力遅延パラメーター nk によって B 多項式の先頭に含まれる固定のゼロの数が多くなる場合に便利です。max(nk-1,0) のラグを nk から 'IODelay' の値に移動して括り出すことができます。

    たとえば、2 つの入力をもつシステムで、1 つ目の入力には 3 サンプルの遅延があり、2 つ目の入力には 6 サイクルの遅延があるとします。また、それらの入力に対する B 多項式の次数は n であると仮定します。これらの遅延は以下を使用して表現できます。

    • nk = [3 6] — この場合、B 多項式は [0 0 0 b11 ... b1n][0 0 0 0 0 0 b21 ... b2n] になります。

    • nk = [3 6]'IODelay',3 — この場合、B 多項式は [b11 ... b1n][0 0 0 b21 ... b2n] になります。

ノイズ チャネルへの積分器の追加。'IntegrateNoise' と長さ Ny の logical ベクトルで構成されるコンマ区切りのペアとして指定します。Ny は出力の数です。

特定の出力に対して 'IntegrateNoise'true に設定すると、そのチャネルに対する ARIX モデルまたは ARI モデルが作成されます。ノイズの積分は、外乱が非定常である場合に便利です。

'IntegrateNoise' を使用する場合、出力チャネルのデータも積分する必要があります。例については、ARIX モデルを参照してください。

出力引数

すべて折りたたむ

推定データを適合する ARX モデル。離散時間の idpoly オブジェクトとして返されます。このモデルは、指定したモデル次数、遅延、および推定オプションを使用して作成されます。

推定結果と使用されたオプションに関する情報は、モデルの Report プロパティに格納されます。Report には次のフィールドがあります。

Report のフィールド説明
Status

モデルのステータスの概要。モデルが構築によって作成されたものか推定によって取得されたものかを示します。

Method

使用された推定コマンド。

InitialCondition

モデル推定時の初期条件の処理。次の値のいずれかとして返されます。

  • 'zero' — 初期条件をゼロに設定。

  • 'estimate' — 初期条件を独立した推定パラメーターとして処理。

このフィールドは、推定オプション セットで InitialCondition オプションが 'auto' に設定されている場合に、初期条件がどのように扱われたかを確認するのに特に便利です。

Fit

推定の定量的評価。構造体として返されます。これらの品質メトリクスの詳細については、Loss Function and Model Quality Metricsを参照してください。構造体には、以下のフィールドがあります。

フィールド説明
FitPercent

正規化平方根平均二乗誤差 (NRMSE)。モデルの応答が推定データにどの程度適合するかをパーセンテージで示す尺度で、fitpercent = 100(1-NRMSE) として表されます。

LossFcn

推定完了時の損失関数の値。

MSE

平均二乗誤差 (MSE)。モデルの応答が推定データにどの程度適合するかを示す尺度です。

FPE

モデルの最終予測誤差。

AIC

生の赤池情報量基準 (AIC)。モデルの品質を示す尺度です。

AICc

小さいサンプルサイズの補正された AIC。

nAIC

正規化された AIC。

BIC

ベイズ情報量基準 (BIC)。

Parameters

モデル パラメーターの推定値。

OptionsUsed

推定に使用されたオプション セット。これは、カスタム オプションを構成していない場合は既定のオプションのセットになります。詳細については、arxOptions を参照してください。

RandState

推定開始時の乱数ストリームの状態。推定時にランダム化が使用されなかった場合は空 [] になります。詳細については、rng を参照してください。

DataUsed

推定に使用されたデータの属性。次のフィールドをもつ構造体として返されます。

フィールド説明
Name

データ セットの名前。

Type

データ型。

Length

データ サンプルの数。

Ts

サンプル時間。

InterSample

入力サンプル間動作。次の値のいずれかとして返されます。

  • 'zoh' — ゼロ次ホールドにより、サンプル間で区分的に一定な入力信号を維持。

  • 'foh' — 1 次ホールドにより、サンプル間で区分的に線形な入力信号を維持。

  • 'bl' — 帯域幅を制限した動作により、ナイキスト周波数を超える連続時間入力信号のパワーがゼロになるように指定。

InputOffset

推定時に時間領域の入力データから削除されたオフセット。非線形モデルの場合は [] になります。

OutputOffset

推定時に時間領域の出力データから削除されたオフセット。非線形モデルの場合は [] になります。

Report の使用の詳細については、Estimation Reportを参照してください。

推定される初期条件。initialCondition オブジェクトまたは initialCondition 値のオブジェクト配列として返されます。

  • 単一実験データ セットの場合、ic は、推定される初期条件 (x0) に対する伝達関数モデルの自由応答 (行列 A および C) を状態空間形式で表します。

  • Ne 回の実験用の複数実験データ セットの場合、ic は、実験ごとに 1 つの initialCondition 値のセットを含む長さ Ne のオブジェクト配列です。

詳細については、initialCondition を参照してください。この引数の使用例については、初期条件の取得を参照してください。

詳細

すべて折りたたむ

ARX 構造

ARX モデルの名前は "追加入力を伴う自己回帰" を意味します。AR モデルとは異なり、ARX モデルには入力項があるため、このように呼ばれています。入力項が外生変数である ARX は "外生変数を伴う自己回帰" とも呼ばれます。ARX モデル構造は次の方程式で与えられます。

y(t)+a1y(t1)+...+anay(tna)=b1u(tnk)+...+bnbu(tnbnk+1)+e(t)

パラメーター na および nb は ARX モデルの次数、nk は遅延です。

  • y(t) — 時間 t における出力

  • na — 極の数

  • nb — ゼロの数

  • nk — 入力が出力に影響する前のシステムの "むだ時間" とも呼ばれる時間に発生する入力サンプルの数

  • y(t1)y(tna) — 現在の出力が依存している前の出力

  • u(tnk)u(tnknb+1) — 現在の出力が依存している前の入力と遅延入力

  • e(t) — ホワイトノイズ外乱の値

よりコンパクトに差分方程式を記述すると次のようになります。

A(q)y(t)=B(q)u(tnk)+e(t)

"q" は遅れ演算子です。具体的には次のようになります。

A(q)=1+a1q1++anaqna

B(q)=b1+b2q1++bnbqnb+1

ARIX モデル

ARIX (追加入力を伴う自己回帰和分) モデルは、ノイズ チャネルに積分器を含む ARX モデルです。ARIX モデル構造は次の方程式で与えられます。

A(q)y(t)=B(q)u(tnk)+11q1e(t)

ここで、11q1 はノイズ チャネル e(t) の積分器です。

AR 時系列モデル

入力がなく、出力が 1 つで、A 多項式の次数が na の時系列データの場合、モデルは次数 na の AR 構造をもちます。

AR (自己回帰) モデル構造は次の方程式で与えられます。

A(q)y(t)=e(t)

ARI モデル

ARI (自己回帰和分) モデルは、ノイズ チャネルに積分器を含む AR モデルです。ARI モデル構造は次の方程式で与えられます。

A(q)y(t)=11q1e(t)

多入力単出力モデル

nu 個の入力をもつ多入力単出力システム (MISO) の場合、nb と nk は行ベクトルです。行ベクトルの i 番目の要素が、列ベクトル u(t) の i 番目の入力に関連付けられている次数と遅延に対応します。同様に、B 多項式の係数も行ベクトルです。ARX MISO 構造は次の方程式で与えられます。

A(q)y(t)=B1(q)u1(tnk1)+B2(q)u2(tnk2)++Bnu(q)unu(tnknu)

多入力多出力モデル

多入力多出力システムの場合、nanb、および nk に出力信号ごとに 1 つの行が格納されます。

多出力のケースでは、arx により、予測誤差の共分散行列 (ノルム) のトレースが最小化されます。

t=1NeT(t)e(t)

重み行列 Lambda を使用して、このノルムを任意の二次ノルムに変換します。

t=1NeT(t)Λ1e(t)

次の構文を使用します。

opt = arxOptions('OutputWeight',inv(lambda))
m = arx(data,orders,opt)

初期条件

時間領域データの場合は、非測定信号が予測子で必要にならないように信号がシフトされます。したがって、初期条件を推定する必要はありません。

周波数領域データの場合、循環畳み込みをサポートする初期条件によるデータの調整が必要になることがあります。

'InitialCondition' 推定オプション (arxOptions を参照) を次の値のいずれかに設定します。

  • 'zero' — 調整なし

  • 'estimate' — 循環畳み込みをサポートする初期条件でデータに対する調整を実行

  • 'auto' — データに基づいて 'zero' または 'estimate' を自動で選択

アルゴリズム

最小二乗の推定問題を構成する一連の過決定の線形方程式については QR 分解で解きます。

正則化を使用しない場合、ARX モデルのパラメーターのベクトル θ は次の正規方程式を解くことで推定されます。

(JTJ)θ=JTy

ここで、J はリグレッサー行列、y は測定出力です。したがって、次のようになります。

θ=(JTJ)1JTy

正則化を使用する場合は正則化項が追加されます。

θ=(JTJ+λR)1JTy

ここで、λ と R は正則化定数です。正則化定数の詳細については、arxOptions を参照してください。

回帰行列が arxOptions で指定された MaxSize よりも大きい場合、データがセグメント化され、データ セグメントに対して反復して QR 分解が実行されます。

バージョン履歴

R2006a より前に導入

すべて展開する